method of lines for coupled fisher kolmogorov pdes

I am trying to figure out how to solve two coupled pdes with method of lines, I did it with pdepe but I want to create my own mesh so i can work with the coupled pdes as stochastic pdes. My pdes are fisher kolmogrov pdes and what i mainly have problems is with the indices and how to break one vector into two so i can get both pdes then add them back together. Here is what i have so far(there might be a lot of it wrong) and there is errors i have.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% u_t= D1*u_xx + gamma1*u(1-u-B1*v)+c1*u_x %%%
%%% v_t= D2*v_xx + gamma2*v(1-v-B2*u)+c2*v_x %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% EXACT Solution %%%%%%%%%%%%%%%%%%%%
%%%%%% MI =(1+m1*exp(n1*(x-w1*t))).^(-2) %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% by finite difference(central) %%%%%%%%%%%
%%% u_t = alpha*u_xx %%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% u_x = (u_{i+1} - u_{i-1}) / 2*h %%%%%%%%%%%%
%%% u_xx = (u_{i+1} - 2u_i + u_{i-1))/h^2) %%%%%
%%% u(i)_t = (u_{i+1}-2u_i+u_{i+1))/h^2 %%
%%%%%%%%%% -gamma1*u(1-u-B1*v) %%%%%%%%%%%%%%%%%
%%%%%%%%%% +c1*((u_{i+1} - u_{i-1}) / 2*h) %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all
clc

global D1 D2 gamma1 gamma2 B1 B2 c1 c2 m1 m2 n1 w1 n2 w2 r n m h
%constant parameters
D1=1;
D2=1;
gamma1=1;
gamma2=1;
B1=1;
B2=2-B1;
c1=1;
c2=2;
m1=1;
m2=1;
n1= -sqrt((1-B1)/6);
w1=-5*sqrt((1-B1)/6)-1;
n2=sqrt((B2-1)/6);
w2=5*sqrt((B2-1)/6)-1;
r=2;

%space discretization
%time and space discretization
n=50;
a=-10;
b=10;
m=50;
t0=0;
T=10; 

vx=linspace(a,b,n);  h=vx(2)-vx(1);
vt=linspace(t0,T,m);   k=vt(2)-vt(1);
%[T,X]=meshgrid(vt,vx);

%Initial condition
u0=(1+m1*exp(n1*vx)).^(-r);
v0=(1+m2*exp(n2*vx)).^(-r);
du0=-(2*exp(n1*vx)*n1)./(1+exp(n1*(vx))).^(r+1);
dv0=-(2*exp(n2*vx)*n2)./(1+exp(n2*(vx))).^(r+1);

%y0=[u0,v0];
y0=[u0,v0,du0,dv0];
%u=y(:,:,1)
%v=y(:,:,2)
%total=u+v


[t,y] = ode45(@(t,y) scfk2(t,y,D1,D2,gamma1,gamma2,c1,c2,B1,B2,h,n,m),vt,y0);
%figure(1)
plot(t,y(:,1:10),'r-o',t,y(:,11:20),'b')
%legend('u','v')
%figure(2)
%surf(vx,vt,y(:,end))
%figure(3)
%surf(vx,vt,y(:,:,2))
%%
function dydt= scfk2(t,y,D1,D2,gamma1,gamma2,c1,c2,B1,B2,h,n,m)

g1=@(y) y;
g2=@(y) 2*y;
ua=@(t) 1;
ub=@(t) 0;
va=@(t) 0;
vb=@(t) 1;

dydt=zeros(n,m);
dudt=dydt(1:n-1); %n-1 entries
dvdt=dydt(n:2*n-2); %n-1 entries

dudt(:,n)   = g1(y);
dudt(1,:)   = ua(t);
dudt(end,:) = ub(t);
dvdt(:,1)   = g2(y);
dvdt(1,:)   = va(t);
dvdt(end,:) = vb(t);


%FTCS solving for dudt and dvdt
for j=1:n-1
    for i=2:m-1
       %Central finite with step size h
        u_x(i,j)= (u(i+1,j)-u(i-1,j))./(2*h);
        u_xx(i,j)=(u(i+1,j)-2*u(i,j)+u(i-1,j))./(h.^2);
        v_x(i,j)=(v(i+1,j)-v(i-1,j))./(2*h);
        v_xx(i,j)=(v(i+1,j)-2*v(i,j)+v(i-1,j))./(h.^2);
        
        %The two pdes
        dudt(i,j) = D1*u_xx(i,j)+gamma1*u(i,j)*(1-u(i,j)-B1*v(i,j))+c1*u_x(i,j);
        dvdt(i,j)= D2*v_xx(i,j)+gamma2*v(i,j)*(1-u(i,j)-B2*v(i,j))+c2*v_x(i,j);
end
    dydt=[dudt;dvdt]; %adding the vectors back to one vector
end
end

Answer 1

[darova]

see this example. It should be helpfull