How to plot a 3D surface from scatter data inside a loop?

Hello there, this is my first post here, and i'd like to thank you all in advance. So my problem is that i'm trying to plot a 3D surface from the scatter data result of an ODE inside a loop. I'm using this to understand how this variables interfere in a electric field used in a quantum control scenario. The thing is that the relation between the variables is not direct, ie, the X and Y components are inside a loop and they are used to create an expression that is interpolated and than used to solve the ODE. I can plot the individual points using Scatter3, but i can't use those points to make a surface. The code i'm using is below. How should i proceed? Thanks again!

    % Variables
        ft = -20000:.1:20000;
        for alpha = 0.001:0.0005:0.01
        for beta = 0.001:0.0005:0.01 
        ai = 0.4;         
        ai2 = 0.6;
        af = 1; 
        w0 = 0.02;
        mi = 6; 
        phii = 0;
        phif = pi/3.2;
    % Functions
        g = 1./(1+exp(-alpha.*ft));
        f = ai*(1-g)+af*g;
        p = 1./(1+exp(-beta.*ft));
        h = phii*(1-p)+phif*p;
    % Electric Field
        E = alpha.*(af-ai).*exp(alpha.*ft).*sin(w0.*ft+h)./(mi.*(1+exp(alpha.*ft)).*sqrt((1-ai+(1-af).*exp(alpha.*ft)).*(ai+af.*exp(alpha.*ft))))+2.*beta.*(phif-phii).*exp(beta.*ft).*sqrt(f.*(1-f)).*cos(w0.*ft+h)/(mi.*(1-2.*f).*(1+exp(beta.*ft).^2));
    % ODE
        [t,y] = ode45(@(t,y) myode(t,y,ft,E,mi,w0), ft, [sqrt(ai)*exp(1i.*0) sqrt(ai2)*exp(1i.*0)]);
        z=1-abs(y(400000,1)).^2;
    % Plot
        scatter3(alpha,beta,z)
        xlabel('Alpha')
        ylabel('Beta')
        hold on

        N = 50;
        xi = linspace(min(alpha),max(alpha),N);
        yi = linspace(min(beta),max(beta),N);
        [X,Y] = meshgrid(xi,yi);
        Z = griddata(alpha,beta,z,X,Y);
        surf(X,Y,Z)
        end
        end
        
% End of the program
function dydt = myode(t,y,ft,E,mi,w0)
E = interpn(ft,E,t);
dydt = [1i.*mi.*y(2).*E.*exp(-1i.*w0.*t);1i.*mi.*y(1).*E.*exp(1i.*w0.*t)];
end 

Answer 1

[Star Strider]

Your data are gridded, however you need to use the reshape function to create matrices from them. I changed your code slightly to save the relevant variables to a matrix (after preallocating ‘abz’ before the loop, and introducing a counter variable ‘kntr’ at the start of the ‘beta’ loop, and took the scatter3 call completely out of the loops):

        for beta = 0.001:0.0005:0.01 
            kntr = kntr + 1

then:

    % ODE
        [t,y] = ode45(@(t,y) myode(t,y,ft,E,mi,w0), ft, [sqrt(ai)*exp(1i.*0) sqrt(ai2)*exp(1i.*0)]);
        z=1-abs(y(400000,1)).^2;
        abz(kntr,:) = [alpha, beta, z];        

and then after the loops:

        ABZ = table(abz(:,1), abz(:,2), abz(:,3), 'VariableNames',{'alpha','beta','z'})
        writetable(ABZ, fullfile(dirpath, 'ABZ_20190408.txt'))

with ‘dirpath’ being the directory you want to write the file to. I attached the file to my Answer.

The plot is then created as:

ABZ = readtable('ABZ_20190408.txt');

alphamtx = reshape(ABZ.alpha, 19, []);
betamtx = reshape(ABZ.beta, 19, []);
zmtx = reshape(ABZ.z, 19, []);

N = 50;
xi = linspace(min(ABZ.alpha),max(ABZ.alpha),N);
yi = linspace(min(ABZ.beta),max(ABZ.beta),N);
[X,Y] = meshgrid(xi,yi);
Z = griddata(alphamtx,betamtx,zmtx,X,Y);

figure
surfc(X,Y,Z)
grid on
xlabel('\alpha')
ylabel('\beta')
zlabel('z')
shading('interp')

And your complete revised code (except for the table and writetable calls) is:

% Variables
kntr = 0;
abz = zeros(19^2, 3);
ft = -20000:.1:20000;
for alpha = 0.001:0.0005:0.01
    for beta = 0.001:0.0005:0.01
        kntr = kntr + 1
        ai = 0.4;
        ai2 = 0.6;
        af = 1;
        w0 = 0.02;
        mi = 6;
        phii = 0;
        phif = pi/3.2;
        % Functions
        g = 1./(1+exp(-alpha.*ft));
        f = ai*(1-g)+af*g;
        p = 1./(1+exp(-beta.*ft));
        h = phii*(1-p)+phif*p;
        % Electric Field
        E = alpha.*(af-ai).*exp(alpha.*ft).*sin(w0.*ft+h)./(mi.*(1+exp(alpha.*ft)).*sqrt((1-ai+(1-af).*exp(alpha.*ft)).*(ai+af.*exp(alpha.*ft))))+2.*beta.*(phif-phii).*exp(beta.*ft).*sqrt(f.*(1-f)).*cos(w0.*ft+h)/(mi.*(1-2.*f).*(1+exp(beta.*ft).^2));
        % ODE
        [t,y] = ode45(@(t,y) myode(t,y,ft,E,mi,w0), ft, [sqrt(ai)*exp(1i.*0) sqrt(ai2)*exp(1i.*0)]);
        z=1-abs(y(400000,1)).^2;
        abz(kntr,:) = [alpha, beta, z];
    end
end

[guilherme stahlberg]

Thank you so much good sir Star Strider, that really helped me a lot, it worked perfectly! Cheers

[Star Strider]

As always, my pleasure!