Problem with the dimension of an integral and plotting it.

Hi, I'm starting to learn matlab, so by advance thanks for reading me.

I want to integrate a function f(q(E),x) respect to "q" between the points q=0 to q=k(E), then evaluate the answer from x=0 to x=56 and finally plot in a surface f(E,x) , the problems that I get are:

1. During the integral procedure, after evaluate 'Int(f(i)q,0,k(i))', Matlab in each iteration do 'Int(f(i-1)q,0,k(i-1))'.
2. The plot using "fplot" is Ok
3. After extract the data from each plot I get a double array but for each energy the array has different length.
4. 'fsurface' code does not work, I think that maybe it is because the previous item.

Hint: I want to reproduce the fig2(c), of this paper: https://www.fkf.mpg.de/52278/kk163.pdf

%Constans
a    = 56 ;               %Step width [A]
hbar = 1.054*10^(-34);    %Plack constant [J*s]
me   = 9.10938*10^(-31);  %electron mass [Kg]
e    = 1.6*10^(-19);      %electron charge [C] 
meff = 0.4*me ;           %effective mass [me]
rl   = 0.35 ;             %LEFT step - coherent reflection amplitude 
rr   = 0.35 ;             %RIGHT step - coherent reflection amplitude 
phil = -pi ;              %LEFT step - coherent reflection phase shift
phir = -pi ;              %RIGHT step - coherent reflection phase shift
eVtoJ = 1.602*10^-19;     % equvalence between eV y J
E0    = -65*10^(-3)*eVtoJ;%Base Energy in "J"

%Variables and length 
E = linspace(-50*10^(-3)*eVtoJ,100*10^(-3)*eVtoJ,5); % 5 values of energy
n = length(E);

x = linspace(0,56,100);% x-values where I want to evaluate the result of the integral 
m = length(x);

k = sqrt(2*meff*(E-E0)./(hbar)^2)*10^-10; % this is the upper limit of the integral

syms q x 

%Function to integrate
f = (1./(sqrt(k.^2-q.^2)*(1+(rl*rr)^2-2*rl*rr*cos(2*q.*a+phil+phir))))*...
    (((1-rl^2)*(1+rr^2+2*rr*cos(2*q.*(x-a)-phir))+(1-rr^2)*(1+rl^2+2*rl*cos(2*q.*x+phil)))); %5 funciones f(q,x)

%Integral code
 for i=1:n
     g(i,:) = int(f(i),q,0,k(i))
 end

% Code to save each figure in this case thear are 5 figures (5 values of Energy) 
 for t=1:n
   h=figure;
    fplot(g(t),[0 56])
    saveas(h,sprintf('FIG%d.fig',t)); % will create FIG1, FIG2,...
end

%code to extract in two columns each pair of data from figures
open FIG1.fig
A=get(gca,'Children');
x1 = get(A,'xdata');
y1 = get(A,'ydata');
 A=[];
 A(:,1)=x1;
 A(:,2)=y1;
 dlmwrite('FIG1.txt',A,'.');

open FIG2.fig
B=get(gca,'Children');
x2 = get(B,'xdata');
y2 = get(B,'ydata');
 B=[];
 B(:,1)=x1;
 B(:,2)=y1;
 dlmwrite('FIG1.txt',B,'.');

open FIG3.fig
C=get(gca,'Children');
x3 = get(C,'xdata');
y3 = get(C,'ydata');
 C=[];
 C(:,1)=x3;
 C(:,2)=y3;
 dlmwrite('FIG1.txt',C,'.');

open FIG4.fig
D=get(gca,'Children');
x4 = get(D,'xdata');
y4 = get(D,'ydata');
 D=[];
 D(:,1)=x4;
 D(:,2)=y4;
 dlmwrite('FIG1.txt',D,'.');

open FIG5.fig
F=get(gca,'Children');
x5 = get(F,'xdata');
y5 = get(F,'ydata');
 F=[];
 F(:,1)=x5;
 F(:,2)=y5;
 dlmwrite('FIG1.txt',F,'.');