The code I've beeing working in:
% Modelo de Sel'kov en 1D
% Este código en MATLAB simulara el modelo de Sel'kov para la glucólisis en 1D
clear all
close all
%%%%%%%%%%%%%%%%%%%%%
%%%Inicialización %%%
%%%%%%%%%%%%%%%%%%%%%
% Parámetros
bc = 0.6;
ac = 0;
Du = 1;
Dv = 2; %Constantes difusivas
% Información inicial y mallado
% w = 10; %sin patrón
w = 80; % patrón
Nx = 10; % Total de puntos discretizados en el dominio [0,L]
x = linspace(0,w,Nx);
dx = x(2) - x(1);
dt = 1; % tamaño del paso
t = 0:dt:50;
Nt = length(t); %numero de puntos en el tiempo
% Condiciones para la superficie
[X, T] = meshgrid(x, t);
U = 0*X;
V = 0*X;
% Vectores columna (más fáciles)
x = x(:);
t = t(:);
%Condición inicial: pequeña perturbación lejos del estado estable
u = bc*ones(length(x),1) + 0.5*rand(Nx, 1);
v = (bc/(ac + bc^2))*ones(length(x),1) + 0.5*rand(Nx,1);
% Condiciones iniciales salvadas.
U(1,:) = u;
V(1,:) = v;
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Generando la matriz%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
% Usando un esquema implícito
% For u
a = (1+2*Du*dt/dx^2); % valores en la diagonal
b = -Du*dt/dx^2; % valores fuera de la diagonal
main = a*sparse(ones(Nx,1));
off = b*sparse(ones(Nx-1,1));
Bu = diag(main) + diag(off,1) + diag(off,-1); %Matriz dispersa
Bv = diag(main) + diag(off,1) + diag(off,-1); %Matriz dispersa
% Condición de frontera flujo cero
Bu(1, end-1) = -b;
Bu(end, 2) = -b;
Bv(1, end-1) = -b;
Bv(end, 2) = -b; %%%%%%%%%%%%%%%%%%%%%%%%%
%%% Resultados %%%
%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
plot(x,u,'g.-', 'linewidth',1);
hold on;
plot(x,v,'r.-', 'linewidth',1);
hold off;
axis([-1 80 -.01 15.01])
for j = 1:Nt-1
% f y g son los términos no lineales en el modelo de Sel'kov
f= dt.*(-u*+ac.*v+v.*u.^2); % f= dx.^2.*(-u*+ac.*v+v.*u.^2);
g= dt.*(bc-ac.*v-v.*u.^2); % g= dx.^2.*(bc-ac.*v-v.*u.^2);
%f = u.^2./v-bc*u;
%g = u.^2 - v; % En cada paso resolvemos el sistema de ecuaciones
u = Bu\(u + f); %u = Bu\(u + dt.*f);
v = Bv\(v + g); %v = Bv\(v + dt.*g); % Gráficas
plot(x,u,'g.-', 'linewidth',1);
hold on;
plot(x,v,'r.-', 'linewidth',1);
hold off;
axis([-1 80 -.01 15.01])
title(['t = ', num2str(j*dt)],'fontsize',24)
drawnow; % Datos para la superficie
U(j,:) = u;
V(j,:) = v;
end %%%% Gráfica de la superficie %%%%
figure(1);
s = surf(x, t, U);
set(s, 'EdgeColor', 'none', 'FaceColor', 'interp');
xlabel('<----x---->')
ylabel('<----t---->')
zlabel('<----u---->') figure(2);
s = surf(x, t, -V);
set(s, 'EdgeColor', 'none', 'FaceColor', 'interp');
xlabel('<----x---->')
ylabel('<----t---->')
zlabel('<----u---->') %%%% contour plot %%%
figure(3);
p = pcolor(x, t, U);
set(p, 'EdgeColor', 'none', 'FaceColor', 'interp');
figure(4);
q = pcolor(x, t, -V);
set(q, 'EdgeColor', 'none', 'FaceColor', 'interp');
The problem is in line 52 because I can't build the matrix in the form I wish.
I need to build those matrices in order to solve $u^{n+1}=A^{-1}(Bu^n+f)$ and $v^{n+1}=A^{-1}(Bu^n+g)$ (line 84 and 85). I can't build the matrices in the form explained in the image.
How can I make the modification in order to obtain what I want?
