Feedback and help with my code that uses ODE15s

Hi all, I've developed some code that solves a moving boundary PDE with ODE15s, that I'd quite like some feedback on as it's still not producing the answers I would expect it to do from the available theory. In order to be able to use ODE15s, I have spatially discretised using the method of lines. On the fixed domain , the Original PDE system is;`

where . For the initial and boundary conditions, we have , , and . I've created two functions, one that solves for S and n using ODE15s called CellsODE, and then one that simulates called CellVelocity. Due to stability conditions, I've implemented an upwind finite difference scheme where the a forward discretisation is used for the derivative in the first equation. The other spatial second order derivatives use a central discretisation, with a forward and backward one at each boundary respectively. My velocity function is,

%CELLVELOCITY Finds V vector given n,s and parameters
function V = CellVelocity(n,s,dx) 
nm=[NaN;n(1:end-1)]; %n_{i-1}
np = [n(2:end);NaN]; %n_{i+1}

dndx=[NaN;np(2:(end-1))-nm(2:(end-1));3*n(end)-4*n(end-1)+n(end-2)]/2*dx; %dndx derivative

V=[0;-dndx(2:end)]/s; %Definition of the velocity, V=0 at the boundary
end

The CellsODE function (for the first and last equation is)

%Extract n and s from the input.
s=input(space_steps+1);
n=[input(1:space_steps-1);0]; %n=0 at the boundary
x=linspace(0,1,space_steps).'; %discrete mesh 

V = CellVelocity(n,s,dx); %Calling for the velocity

np = [n(2:end);NaN]; %n_{i+1} %Declaring the indexes
vm = [NaN;V(1:end-1)]; %v_{i-1}
vp = [V(2:end);NaN]; %v_{i+1}

dndx = [np(1:end-1)-n(1:end-1);n(end)-n(end-1)]/dx; %dndx derivative, first order forward 
dvdx = [-V(3)+4*V(2)-3*V(1);vp(2:end-1)-vm(2:end-1);3*V(end)-4*V(end-1)+V(end-2)]/2*dx; %dvdx derivative, second order

dndt = (x*V(end)-V).*dndx/s   -n.*dvdx/s  +n.*(1-n);  %ode for n 
dsdt = V(end); %ode for s

output=[dndt;dsdt];
end

The 'Command script' that defines all the variables and calls CellsODE etc is;

T = 1000000 ; %Time
space_steps = 200; %Space steps
x = linspace(0,1,space_steps).'; %x vector
dx = 1/(space_steps-1); %space step
s0 = 1; %initial s
n0 = 0.6*(1-x.^2); %initial conditions
T_span=[0;T]; %T_vector

[T_out,Y_out] = ode15s(@(t,y) CellsODE(t,y,space_steps,dx),T_span,[n0;s0]);

My main interest is finding the function , and one can plot the adaptvie time mesh against the found by doing

Boundary = Y_out(:,end);
plot(T_out, Boundary);

Theory tells us that S for large times should be linear and the gradient should be equal to , however mine is coming out to be something very small and completely dissimilar. This is the reason why I think there is something wrong with my code, as it does not match up with the theory. One can calculate the gradient there by simply doing

grad = gradient(Boundary,T_out);

I've tried many things, I've tried changing all the discretisations to be other ones, Changing the positioning of the variables S in the code and everything under the sun. What I have noticed is that changing the number of mesh points really changes the output, which I don't think should happen. I've attached an image of the Boundary, one can see it almost goes linear but if you closer then it isn't. Thanks for reading.