MATLAB: Please help me visualize the analytical solution of the model

3d plotsMATLABMATLAB and Simulink Student Suitepdevisualize

I am currently trying to plot this solution to get the visualisation. This is the code I have currenlty and I am getting several errors, and I am not sure this is correct. Can someone help with this problem? 
clc;clear all;
m=1;m0=0;
n=1;n0=0;
t=0;
J = 10;                     
L = 10;
dx = (m-m0) / J;                % dx: mesh size

dy = (n-n0) / L;                % dx: mesh size
tf = 0.1;                   % final simulation time
Nt = 50;                    % Nt: number of time steps
dt = tf/Nt;
k = 0;
v=1;
Ds=1;
Kt=1;
Ci=0;
Ca=0;
Cb=0;
Cc=0;
Cd=0;
alpha=v/(2*Ds);
beta=v/(2*Ds);
gamma=-(v^2/(2*Ds)+Kt);
% Evaluate the initial conditions
x = m0 : dx : m;              % generate the grid point
y = n0 : dy : n; 
% store the solution at all grid points for all time steps
% u = zeros(J+1,Nt);
 u_ex = zeros(J+1,L+1,Nt) ;
% Find the approximate solution at each time step
for n = 1:Nt-1
    t = n*dt;         % current time
      
    % calculate the analytic solution
    
    for i=2:J+1
        xj = m0 + (i-1)*dx;
        for j=2:L+1
            yj=n0+(j-1)*dy;
          
            for p=1:10
                c_ex1(xj,yj)=(sin(p*pi*xj/m)*(((-2)*Cc*p*pi*exp((-1)*gamma*t)*(1-exp((-1)*alpha*m)*cos(p*pi)))/(sinh(p*pi*n/m)*((alpha*m)^2+(p*pi)^2)))*sinh(p*pi*(yj-n)/m)+(((-2)*Cd*p*pi*exp((-1)*beta*n+(-1)*gamma*t)*(1-exp((-1)*alpha*m)*cos(p*pi)))/(sinh(p*pi*n/m)*((alpha*m)^2+(p*pi)^2)))*sinh(p*pi*yj/m))+(sin(p*pi*yj/n)*(((-2)*Cb*p*pi*exp((-1)*alpha*m+(-1)*gamma*t)*(1-exp((-1)*beta*n)*cos(p*pi)))/(sinh(p*pi*m/n)*((beta*n)^2+(p*pi)^2)))*sinh(p*pi*(xj-m)/n)+(((-2)*Ca*p*pi*exp((-1)*gamma*t)*(1-exp((-1)*beta*n)*cos(p*pi)))/(sinh(p*pi*m/n)*((beta*n)^2+(p*pi)^2)))*sinh(p*pi*xj/n));
            
            for s=1:10
                for r=1:10
                   c_ex2(s,r,t)=((4/(m*n))*int(int(-c_ex1(xj,yj)*sin(r*pi*xj/m)*sin(s*pi*yj/n))))*sin(r*pi*xj/m)*sin(s*pi*yj/n)*exp(-Ds((r^2)/m+(s^2)/n)*pi^2*t); 
                end
            end
            end
            u_ex=exp(alpha*xj+beta*yj+gamma*t)*(c_ex1(p)+c_ex2(s,r,t));
        end
    
    end
    
    
end
% Plot the results
tt = dt : dt : Nt*dt;
figure(1)
surf(x,y,tt, u_ex');     % 3-D surface plot
xlabel('x')
ylabel('t')
zlabel('u')
title('Analytic solution of 2-D parabolic equation')
i am trying to figure out this solution
where
and

Best Answer

Thank you for clear code. Everything looks correct (i didn't tested. it only looks)
You have x y and t independent variables, it means that you should have 3D matrix
u_ex(i,j,n) = exp(alpha*xj+beta*yj+gamma*t)*(phixy+PHIxyt);
This part with double integral
for s=1:10
    for r=1:10
       r1=Ci*exp(-alpha*xj-beta*yj)-phixy
       r2=sin(r*pi*xj/m)
       r3=sin(s*pi*yj/n)
       r4=exp(-Ds*((r^2)/m+(s^2)/n)*pi^2*t)      
       PHIxyt=((4/(m*n))*int(int( r1*r2*r3,xj,0,m),yj,0,n))*r2*r3*r4; 
    end
end
It needed to be re-written (integral2)
for s=1:10
    for r=1:10
       r1 = @(xj,yj) Ci*exp(-alpha*xj-beta*yj)-phixy
       r2 = @(xj) sin(r*pi*xj/m)
       r3 = @(yj) sin(s*pi*yj/n)
       r4 = exp(-Ds*((r^2)/m+(s^2)/n)*pi^2*t)      
       R = @(x,y) r1(x,y)*r2(x)*r3(x);
       PHIxyt = ((4/(m*n))*integral2( R,0,m,0,n ) *r2(xj)*r3(yj)*r4; 
    end
end
To visualize all this i'd use isosurface
isosurface(u_ex,val)
I didn't tested all those things on my computer (i'm just scared). Let me know if there is something wrong
Related Question