Schnakenberg system

clear all ;
 close all;


D1=1;
D2=.1;

N=200; 
L=100; 
c=0;

t_end=100;

x =linspace (c,L,N);
dx= x(2) -x(1);
dt=.001         
t =1:dt:t_end;

u = zeros(length(x),length(t));
v = zeros(length(x),length(t));


 u(:,1) = 0;

 u(70:140,1) = 150;
 v(:,1) =0;

 v(70:140,1) = 150;

for n=1:length(t)-1
   
    for i=2:length(x)-1
        u(i,n+1) = u(i,n) + D1*dt/(dx^2)*(u(i-1,n)-2*u(i,n)+u(i+1,n));%+dt*a-dt*T_mat1(i,n)+dt*(T_mat1(i,n)^2*(T_mat2(i,n)));
        v(i,n+1) = v(i,n) + D2*dt/(dx^2)*(v(i-1,n)-2*v(i,n)+v(i+1,n));%+dt*b-dt*(T_mat1(i,n)^2*(T_mat2(i,n)));
    end
    
    u(1,n+1) = u(1,n) +dt/2*(dx)*(u(2,n)-u(1,n));    
    v(1,n+1) = v(1,n) +dt/2*(dx)*(v(2,n)-v(1,n)); 
    
     u(end,n+1)= u(end,n)+dt/2*(dx)*(u(end-1,n)-u(end,n));
     v(end,n+1)= v(end,n)+dt/2*(dx)*(v(end-1,n)-u(end,n));

end


  for n=1 :length(t)
    Total_u(n)= sum(u(1:length(x)-1,n))*dx;
  end
 
  for n=1 :length(t)
    Total_v(n)= sum(v(1:length(x)-1,n))*dx;
 end
   

figure
[tt,xx] = meshgrid (t,x);
zz=tt.*xx;
mesh(xx,tt,u);
xlabel('X coordinate (m)')
ylabel('Time (sec)')
zlabel('Temperature (F)')
zlim([0,200])

figure
[tt,xx] = meshgrid (t,x);
zz=tt.*xx;
mesh(xx,tt,v);
xlabel('X coordinate (m)')
ylabel('Time (sec)')
zlabel('Temperature (F)')
zlim([0,200])

figure

hold on
plot(x, u(:,1),x,v(:,1)) % initial condition T(x,0) vs x
% plot(x_vec, T_mat(:, 500)) % T(x, t_end/3) vs x
%plot(x_vec, T_mat(:, round(length(t_vec)/4)))
plot(x, u(:, round(length(t)/3)),x,v(:, round(length(t)/3))) % T(x, t_end/3) vs x
plot(x, u(:, round(length(t)*(2/3))),x, v(:, round(length(t)*(2/3)))) % T(x, t_end*2/3) vs x
%plot(x_vec, T_mat(:, round(length(t_vec)*(3/4))))
plot(x, u(:, length(t)),x,v(:, length(t))) % final state T(x, t_end) vs x
legend('t=0','t=(1/3)t_{end}', 't=(2/3)t_{end}','t=t_{end}')
hold off