# Applying the method of lines to a partial differential equation and using Runge-kutta method in a matlab code

By method of lines I converted the PDE

u_t=u_{xx}, with the initial and boundary conditions

u(0,t)=u(0.1,t)
u(0,x)=sin(2*pi*x)


to a system of ODEs, as follows

function ut=pde1(t,u)
%
% Problem parameters
% global ncall
xl=0.0;
xr=1.0;
d=10;
a=10;
%
% PDE
n=length(u);
h=((xr-xl)/(n-1));
for i=1:n
if(i==1)     ut(i)=0.0;
elseif(i==n) ut(i)=0;
else         ut(i)=d*(u(i+1)-2.0*u(i)+u(i-1))/h^2;
end
end
ut=ut';
%


Now I am solving the ODEs ut by Runge-kutta method as follows:

%Initial condition
n=500;
for i=1:n
u0(i)=sin((2*pi)*(i-1)/(n-1));
end
t0=0.0;
tf=0.1;
tout=linspace(t0,tf,n);
y=zeros(n,length(tout));
y(:,1)=u0';
h=1/n;

for i = 1 : length(tout)
t=tout(i);
k1 = pde_1(t,y(:,i));
k2 = pde_1(t+h/2, y(:,i)+h*k1/2);
k3 = pde_1(t+h/2, y(:,i)+h*k2/2);
k4 = pde_1(t+h, y(:,i)+h*k3);
y(:,i+1) = y(:,i) + h*(k1 + 2*k2 + 2*k3 + k4)/6;

end

y
plot(tout,y(:,1))


but the outputs are NAN! any help?

• Why don't you use one of the Matlab built-in version of Runge-Kutta algorithm like ode45 ? – Jean Marie Jun 8 '17 at 2:24
• @JeanMarie I am avoiding the build functions at this level for some a reason, I am considering some different types of runge kutta – Raed Jun 8 '17 at 4:36
• Please correct the boundary conditions, order of variables, 0.1 vs. 1.0, is this a periodic b.c. or are both boundaries zero? If periodic, then your implementation is wrong. – LutzL Jun 8 '17 at 10:20
• There are several problems with the code. As @LutzL pointed out, the problem statement is not consistent with the implementation. However, the reason behind your NaNs is that the CFL condition is not satisfied. For the diffusion equation you need the the time step $h$ to scale as $dx^2$ (where $dx$ is the spatial step), but in the code its simply set to $h = dx/10$. – ekkilop Jun 20 '17 at 15:50
• @ekkilop you right, the stability condition should be satisfied – Raed Jun 23 '17 at 18:22