How to incorporate v and w

%% u_t = D1 * u_xx,  u(0,t) = u(1,t) = 0, u(x,0)= cos(pi*(x-0.5)), D1 = 0.1.
% % v_t = D2 * v_xx + a*w ,  v(0,t) = v(1,t) = 0, v(x,0)= sin(pi*(x-0.5)),D2 = 0.2, a = 1.5.
% % w_t = D3 * w_xx + c*u + b*v,  w(0,t) = w(1,t) = 0, w(x,0)= e^(2*x),D3 = 0.3, b = 2, c = 3.
%subdivisions space /time
ht = 0.01; Tmax = 1.2; nx = 33; hx = 1/(nx-1);  D1 = 0.1;  x = [0:hx:1]';
%matrices
K = stiff(D1,hx,nx); M = mass(1/ht,hx,nx);A = M + K;
% disp(A(1))
%boundary conditions by deleting rows
A(nx,:) = [];   A(:,nx) = [];   A(1,:) = [];    A(:,1) = [];
%creation
R = chol(A);
%initial condition
u = cos(pi*(x - 0.5 ));
%time step loop
k = 0;
while (k*ht < Tmax)
    k = k+1;
    %compute the right-hand side
    b = M*u;    b(nx) = [];     b(1) = [];
    %solve
    u = zeros(nx,1);    u(2:nx-1) = R\(R'\b);
    figure(11)
    plot(x,u(:,1),'Linewidth',1.5)
    xlabel \bfx
    ylabel \bfu
    % % % % hold on
end
function R = stiff(nu,h,n)
[m1,m2] = size(nu);
if (length(nu) == 1)
    ee = nu*ones(n-1,1)/h;  e1 = [ee;0]; e2 = [0;ee]; e = e1+e2;
else
    ee = 0.5*(nu(1:n-1)+nu(2:n))/h;  e1 = [ee;0]; e2 = [0;ee]; e = e1+e2;
end
R = spdiags([-e1 e -e2],-1:1,n,n);
return
end
function M = mass(alpha,h,n)
[m1,m2] = size(alpha);
if(length(alpha) == 1)
    ee = alpha*h*ones(n-1,1)/6;
    e1 = [ee;0]; e2 = [0;ee]; e = e1+e2;
else
    ee = h*(alpha(1:n-1) + alpha(2:n))/12;
    e1 = [ee;0]; e2 = [0;ee]; e = e1+e2;
end
M = spdiags([e1 2*e e2],-1:1,n,n);
return
end