# MATLAB: Surface plot of PDE numeric solution

3d plotsMATLABpde

Hello, I'm solving the system of 2 PDE's, each function depends on 3 variables(radius, angle and time). I'm using explicit difference scheme. As a result I've got two 3D matrices (one for each function). To visualise the results, I want to plot surface plots for each function with fixed t(time). For example: In a moment t=0.5 the surface for u(:,:,0.5): X-axis will be the radius, Y-axis will be the angle and Z-axis will be the function value in in this point (x,y). Will be glad for any advices. Thanks.
Here is the code, if it's difficult to understand my English.
if true  % code   clc;clear all;%grid for rn=100;r_min=0.01;r_max=1;hr=(r_max-r_min)/(n-1)r=r_min:hr:r_maxnr=max(size(r))%grid for thetam=100;th_min=0;th_max=2*pi;hth=(th_max-th_min)/(m-1)theta=th_min:hth:th_maxnth=max(size(theta))%grid for tl=10000;t_min=0;t_max=1;ht=(t_max-t_min)/(l-1)time=t_min:ht:t_maxnt=max(size(time))u = zeros(nr,nth,nt);v = zeros(nr,nth,nt);%Inititalizationu0=0;for i=1:nr  for j=1:nth      if (r(i)<r0)          u0=0;      elseif ((r(i)>=r0) &&(r(i)<r1))          u0=(r(i)-r0)/(r1-r0);      elseif ((r(i)>=r1)&&(r(i)<r2))          u0=1;      elseif ((r(i)>=r2)&&(r(i)<r3))          u0=(r(i)-r3)/(r2-r3);      else          u0=0;      end      u(i,j,1)=u0;  endendu(:,:,1)v0=0;for i=1:nr  for j=1:nth      if(r(i)<r1)          v0=1;      elseif ((r(i)>=r1)&&(r(i)<r3))          v0=(r(i)-r3)/(r1-r3);      else          v0=0;      end      v(i,j,1)=v0;  endendv(:,:,1)for t=1:nt-1    for i=2:nr-1        for j=2:nth-1            %derivatives            d2udr2= (u(i+1,j,t)-2*u(i,j,t)+u(i-1,j,t))/hr^2;            dudr=(u(i+1,j,t)-u(i-1,j,t))/(2*hr);            d2udth2=(u(i,j+1,t)-2*u(i,j,t)+u(i,j-1,t))/hth^2;            d2vdr2= (v(i+1,j,t)-2*v(i,j,t)+v(i-1,j,t))/hr^2;            dvdr=(v(i+1,j,t)-v(i,j,t))/hr;            d2vdth2=(v(i,j+1,t)-2*v(i,j,t)+v(i,j-1,t))/hth^2;            %centre nodes            u(i,j,t+1)=u(i,j,t)+ht*(d2udr2+(1/r(i))*dudr+(1/r(i)^2)*d2udth2);%+u(i,j,t)*(1-u(i,j,t)-c*v(i,j,t)));            v(i,j,t+1)=v(i,j,t)+ht*mu*(d2vdr2+(1/r(i))*dvdr+(1/r(i)^2)*d2vdth2);%+a*v(i,j,t)*(1-c*u(i,j,t)-b*v(i,j,t)));        end        %derivatives for left boundary        d2udr2= (u(i+1,1,t)-2*u(i,1,t)+u(i-1,1,t))/hr^2;        dudr=(u(i+1,1,t)-u(i-1,1,t))/(2*hr);        d2udth2=(u(i,2,t)-2*u(i,1,t)+u(i,nth,t))/hth^2;        d2vdr2= (v(i+1,1,t)-2*v(i,1,t)+v(i-1,1,t))/hr^2;        dvdr=(v(i+1,1,t)-v(i-1,1,t))/(2*hr);        d2vdth2=(v(i,2,t)-2*v(i,1,t)+v(i,nth,t))/hth^2;                uleft=u(i,1,t)+ht*(d2udr2+(1/r(i))*dudr+(1/r(i)^2)*d2udth2);%+u(i,1,t)*(1-u(i,1,t)-c*v(i,1,t)));        vleft=v(i,1,t)+ht*mu*(d2vdr2+(1/r(i))*dvdr+(1/r(i)^2)*d2vdth2);%+a*v(i,1,t)*(1-c*u(i,1,t)-b*v(i,1,t)));                %derivatives for right boundary        d2udr2= (u(i+1,nth,t)-2*u(i,nth,t)+u(i-1,nth,t))/hr^2;        dudr=(u(i+1,nth,t)-u(i-1,nth,t))/(2*hr);        d2udth2=(u(i,1,t)-2*u(i,nth,t)+u(i,nth-1,t))/hth^2;        d2vdr2= (v(i+1,nth,t)-2*v(i,nth,t)+v(i-1,nth,t))/hr^2;        dvdr=(v(i+1,nth,t)-v(i-1,nth,t))/(2*hr);        d2vdth2=(v(i,1,t)-2*v(i,nth,t)+v(i,nth-1,t))/hth^2;                uright=u(i,nth,t)+ht*(d2udr2+(1/r(i))*dudr+(1/r(i)^2)*d2udth2);%+u(i,nth,t)*(1-u(i,nth,t)-c*v(i,nth,t)));        vright=v(i,nth,t)+ht*mu*(d2vdr2+(1/r(i))*dvdr+(1/r(i)^2)*d2vdth2);%+a*v(i,nth,t)*(1-c*u(i,nth,t)-b*v(i,nth,t)));         %filling boundaries for theta        u(i,1,t+1)=uleft;        v(i,1,t+1)=vleft;        u(i,nth,t+1)=uright;        v(i,nth,t+1)=vright;    end    %boundaries for r    for j=1:nth        u(1,j,t+1)=u(2,j,t+1);        u(nr,j,t+1)=u(nr-1,j,t+1);        v(1,j,t+1)=v(2,j,t+1);        v(nr,j,t+1)=v(nr-1,j,t+1);    endend  end
surf(r,theta,?)

 clc;clear all;n=100;r_min=0.01;r0=r_min;r_max=1;r3=r_max;r1=r0+0.2;r2=r1+0.2;hr=(r_max-r_min)/(n-1);r=r_min:hr:r_max;nr=max(size(r));%grid for thetam=100;th_min=0;th_max=2*pi;hth=(th_max-th_min)/(m-1);theta=th_min:hth:th_max;nth=max(size(theta));%grid for tl=100;t_min=0;t_max=1;ht=(t_max-t_min)/(l-1);time=t_min:ht:t_max;nt=max(size(time));u = zeros(nr,nth,nt);v = zeros(nr,nth,nt);%Inititalizationu0=0;for i=1:nr  for j=1:nth      if (r(i)<r0 || r(i)>r3)          u0=0;      elseif ((r(i)>=r0) &&(r(i)<r1))          u0=(r(i)-r0)/(r1-r0);      elseif ((r(i)>=r1)&&(r(i)<r2))          u0=1;      elseif ((r(i)>=r2)&&(r(i)<r3))          u0=(r(i)-r3)/(r2-r3);      end      u(i,j,1)=u0;  endendu(:,:,1);v0=0;for i=1:nr  for j=1:nth      if(r(i)<r1)          v0=1;      elseif ((r(i)>=r1)&&(r(i)<r3))          v0=(r(i)-r3)/(r1-r3);      else          v0=0;      end      v(i,j,1)=v0;  endendv(:,:,1);Mu=1;for t=1:nt-1    for i=2:nr-1        for j=2:nth-1            %derivatives            d2udr2= (u(i+1,j,t)-2*u(i,j,t)+u(i-1,j,t))/(hr^2);            dudr=(u(i+1,j,t)-u(i-1,j,t))/(2*hr);            d2udth2=(u(i,j+1,t)-2*u(i,j,t)+u(i,j-1,t))/(hth^2);            d2vdr2= (v(i+1,j,t)-2*v(i,j,t)+v(i-1,j,t))/(hr^2);            dvdr=(v(i+1,j,t)-v(i,j,t))/hr;            d2vdth2=(v(i,j+1,t)-2*v(i,j,t)+v(i,j-1,t))/(hth^2);            %centre nodes            u(i,j,t+1)=u(i,j,t)+ht*(d2udr2+(1/r(i))*dudr+(1/r(i)^2)*d2udth2);%+u(i,j,t)*(1-u(i,j,t)-c*v(i,j,t)));            v(i,j,t+1)=v(i,j,t)+ht*Mu*(d2vdr2+(1/r(i))*dvdr+(1/r(i)^2)*d2vdth2);%+a*v(i,j,t)*(1-c*u(i,j,t)-b*v(i,j,t)));        end        %derivatives for left boundary        d2udr2= (u(i+1,1,t)-2*u(i,1,t)+u(i-1,1,t))/hr^2;        dudr=(u(i+1,1,t)-u(i-1,1,t))/(2*hr);        d2udth2=(u(i,2,t)-2*u(i,1,t)+u(i,nth,t))/hth^2;        d2vdr2= (v(i+1,1,t)-2*v(i,1,t)+v(i-1,1,t))/hr^2;        dvdr=(v(i+1,1,t)-v(i-1,1,t))/(2*hr);        d2vdth2=(v(i,2,t)-2*v(i,1,t)+v(i,nth,t))/hth^2;                uleft=u(i,1,t)+ht*(d2udr2+(1/r(i))*dudr+(1/r(i)^2)*d2udth2);%+u(i,1,t)*(1-u(i,1,t)-c*v(i,1,t)));        vleft=v(i,1,t)+ht*Mu*(d2vdr2+(1/r(i))*dvdr+(1/r(i)^2)*d2vdth2);%+a*v(i,1,t)*(1-c*u(i,1,t)-b*v(i,1,t)));                %derivatives for right boundary        d2udr2= (u(i+1,nth,t)-2*u(i,nth,t)+u(i-1,nth,t))/hr^2;        dudr=(u(i+1,nth,t)-u(i-1,nth,t))/(2*hr);        d2udth2=(u(i,1,t)-2*u(i,nth,t)+u(i,nth-1,t))/hth^2;        d2vdr2= (v(i+1,nth,t)-2*v(i,nth,t)+v(i-1,nth,t))/hr^2;        dvdr=(v(i+1,nth,t)-v(i-1,nth,t))/(2*hr);        d2vdth2=(v(i,1,t)-2*v(i,nth,t)+v(i,nth-1,t))/hth^2;                uright=u(i,nth,t)+ht*(d2udr2+(1/r(i))*dudr+(1/r(i)^2)*d2udth2);%+u(i,nth,t)*(1-u(i,nth,t)-c*v(i,nth,t)));        vright=v(i,nth,t)+ht*Mu*(d2vdr2+(1/r(i))*dvdr+(1/r(i)^2)*d2vdth2);%+a*v(i,nth,t)*(1-c*u(i,nth,t)-b*v(i,nth,t)));         %filling boundaries for theta        u(i,1,t+1)=uleft;        v(i,1,t+1)=vleft;        u(i,nth,t+1)=uright;        v(i,nth,t+1)=vright;    end    %boundaries for r    for j=1:nth        u(1,j,t+1)=u(2,j,t+1);        u(nr,j,t+1)=u(nr-1,j,t+1);        v(1,j,t+1)=v(2,j,t+1);        v(nr,j,t+1)=v(nr-1,j,t+1);    endend figure, surf(r,theta,u(:,:,3)), shading interp, title(' U_{3}'); figure, surf(r,theta,v(:,:,3)), shading interp, title(' V_{3}');