rk4

 RK4 One step of  Explicit Runge-Kutta fourth order scheme.
 
 [u1]=rk4(tt,deltat,f,u0,dx,A,w,visc,uex,uex1,bc);

 Input : tt =time
         deltat = time step
         f =  function for the r.h.s of the parabolic equation
         u0 = initial data
         dx = LGL first derivative matrix
         A = matrix related to differential operator
         w = LGL weights array
         visc = viscosity coefficient (constant > 0)
         uex  = exact solution (uex=@(x)[uex(x)], with .*, .^, ./)
         uex1 = exact solution (uex=@(x)[uex(x)], with .*, .^, ./)
         bc     = choice of boundary conditions: 1 == Dirichlet 
                                                 2 == Neumann 

 Output: u1 = array solution, computed with one step of RK4

0001 function [u1]=rk4(tt,deltat,f,u0,dx,A,w,visc,uex,uex1,cb);
0002 % RK4 One step of  Explicit Runge-Kutta fourth order scheme.
0003 %
0004 % [u1]=rk4(tt,deltat,f,u0,dx,A,w,visc,uex,uex1,bc);
0005 %
0006 % Input : tt =time
0007 %         deltat = time step
0008 %         f =  function for the r.h.s of the parabolic equation
0009 %         u0 = initial data
0010 %         dx = LGL first derivative matrix
0011 %         A = matrix related to differential operator
0012 %         w = LGL weights array
0013 %         visc = viscosity coefficient (constant > 0)
0014 %         uex  = exact solution (uex=@(x)[uex(x)], with .*, .^, ./)
0015 %         uex1 = exact solution (uex=@(x)[uex(x)], with .*, .^, ./)
0016 %         bc     = choice of boundary conditions: 1 == Dirichlet
0017 %                                                 2 == Neumann
0018 %
0019 % Output: u1 = array solution, computed with one step of RK4
0020 %
0021 
0022 %   Written by Paola Gervasio
0023 %   $Date: 2007/04/01$
0024 %
0025 
0026 np=length(u0);
0027 del2=deltat*0.5;
0028 K1=feval(f,tt,deltat,u0,A,dx,w,visc,uex,uex1,cb);
0029 t1=tt+del2; u1=u0+del2*K1;
0030 if cb ==1
0031 u1(1)=uex(-1,t1); u1(np)=uex(1,t1);
0032 end
0033 
0034 K2=feval(f,t1,deltat,u1,A,dx,w,visc,uex,uex1,cb);
0035 u1=u0+del2*K2;
0036 if cb ==1
0037 u1(1)=uex(-1,t1); u1(np)=uex(1,t1);
0038 end
0039 K3=feval(f,t1,deltat,u1,A,dx,w,visc,uex,uex1,cb);
0040 t1=tt+deltat; u1=u0+deltat*K3;
0041 if cb ==1
0042 u1(1)=uex(-1,t1); u1(np)=uex(1,t1);
0043 end
0044 K4=feval(f,t1,deltat,u1,A,dx,w,visc,uex,uex1,cb);
0045 u1=u0+deltat/6*(K1+2*K2+2*K3+K4);
0046 if cb ==1
0047 u1(1)=uex(-1,t1); u1(np)=uex(1,t1);
0048 end
0049 
0050 return
0051