clc
close all
clear all

%Program to simulate the Brusselator

%parameters in a structure
struct.k1A=0.1; %k1*[A]
struct.k2=0.1;
struct.k3B=0.5; %k3*[B]
struct.k4=0.1;
struct.Dx=0; %diffusion for [X]
struct.Dy=0; %diffusion for [Y]

%number of cells
struct.N=1;

%set tspan
tStart=0;
tEnd=1000;

%set minimum and maximum values for initial concentration of [X]
%(initial [X] uniformly distributed random between minX and maxX)
%NOTE: to solve C6 set minX=maxX=desiderated value for X
minX=1;
maxX=1;

%set minimum and maximum values for initial concentration of [Y]
%(initial [Y] uniformly distributed random between minY and maxY)
%NOTE: to solve C6 set minY=maxY=desiderated value for Y
minY=1;
maxY=1;

%generate initial [X] concentration in the N cells
initX=minX+(maxX-minX)*rand(1,struct.N);

%generate initial [Y] concentration in the N cells
initY=minY+(maxY-minY)*rand(1,struct.N);

%call ode45 to solve the system of differential equations
%set the options for ode45
%T: time
%Z: first N elements are time profiles for concentration of [X] in each
%cells; last N elements are time profiles for concentration of [Y] in each
%cells
optionsODE=odeset('AbsTol',1e-8.*ones(1,struct.N*2),'RelTol',1e-8);
[T Z]=ode45('brusselator',[tStart tEnd],[initX initY],optionsODE,struct);

if struct.N==1
    %plot [X] and [Y] in one cell
    X=Z(:,1);
    Y=Z(:,2);
    figure(1)
    plot(T,X,'-b',T,Y,'-r');
    xlabel('time');
    ylabel('concentration');
    legend('[X]','[Y]');
else
    %plot [X] and [Y] in N cells (row: time; column: cell number)
    X=Z(:,1:struct.N);
    Y=Z(:,struct.N+1:struct.N*2);
    figure(1) %3D 
    mesh(T,(1:1:struct.N),X')
    xlabel('time')
    ylabel('cell number')
    zlabel('[X] concentration')
    figure(2)
    mesh(T,(1:1:struct.N),Y')
    xlabel('time')
    ylabel('cell number')
    zlabel('[Y] concentration')
    figure(3) %heat map
    subplot1=subplot(2,1,1);
    imagesc(T,(1:1:struct.N),X')
    title('[X] concentration')
    xlabel('time')
    ylabel('cell number')
    colorbar('peer',subplot1);
    subplot2=subplot(2,1,2);
    imagesc(T,(1:1:struct.N),Y')
    title('[Y] concentration')
    xlabel('time')
    ylabel('cell number')
    colorbar('peer',subplot2);
end