ODE15s doesn't seem to be working, one of the equations seem to be unchanged

I want to solve the following set of equations:

via the finite volume method, so I integrate the above equations over an inteval [h_{i},h_{i+1}] to get a set of ODEs which I then proceed to solve using ODE15s, as this appears to be a stiff set of equations. The attached code seems to be getting the change in u, which is good, BUT, gradients in u should give rise to a RHS for the first equation and therefore evolves nu, but it doesn't. I'm at a loss why. I've checked the fluxes, which should be the source of the issue but they seem fine.

Can anyone point out any trivial errors?

%This is a code to solve the isothermal sintering equations.

S=struct;
%%---Physical parameters---
S.N = 251; %This is the number of cells-1 
S.L = 1; %This is the initial length of the rod
S.g = 9.81;  %Acceleration due to gravity.
S.K = 1e-2; %This is the Laplace constant from the sintering potential.
S.eta_0 = 0.05; %The shear stress of the skeleton
S.T = 8; %This is the time of sintering.
S.rho_m = 1; %This is the density of the individual metallic particles.



%%---Initial conditions
%Length
S.X=linspace(0,S.L,S.N); %This is the partition of the initial length
X = S.X;
%Density
rho_0_int=0.5-0.1*exp(-500*(X-0.5).^2);
rho_0=0.5*(rho_0_int(1:S.N-1)+rho_0_int(2:S.N));
%Amount of mass as a function of position
h_interface=cumtrapz(X,rho_0_int);
S.dh=diff(h_interface);

%Mass at the midpoint;
S.h=0.5*(h_interface(1:S.N-1)+h_interface(2:S.N));

%Total mass
S.M = h_interface(end);

U=S.M/(S.rho_m*S.T); %Scaling for velocity
a_1=S.T/(U*S.M);
a_2=S.eta_0*S.T*S.rho_m/S.M^2;
a_3=S.g*S.T/U;

%Initial velocity
u0=zeros(S.N-1,1);
% Define the system of equations as a function

y0 = [1./rho_0'; u0];
% Finite Volume Method solution loop
tspan = [0 1];  % You can adjust this based on the time range you're interested in
[t, y] = ode15s(@(t,y)(ode_system(t,y,S)), tspan, y0);
rho=1./y(:,1:S.N-1);
u=y(:,S.N:2*S.N-2);
Tspan=length(y(:,1));
disp('PDEs solved, now computing new grid');

PDEs solved, now computing new grid

figure()
hold on
for i=1:Tspan
    plot(u(i,:))
end
hold off

function dydt = ode_system(t, y, S)

U=S.M/(S.rho_m*S.T); %Scaling for velocity
a_1=S.T/(U*S.M);
a_2=S.T*S.rho_m/S.M^2;
a_3=S.g*S.T/U;
    nu = y(1:S.N-1)';  % specific volume
    u = y(S.N:2*(S.N-1))';  % velocity
    %Compute the left and right specific volumes:
    nu_left=nu(1); %15*nu(1)/8-5*nu(2)/4+3*nu(3)/8;
    nu_right=nu(S.N-1); %15*nu(S.N-3)/8-5*nu(S.N-2)/5+3*nu(S.N-1)/8;
    %Compute the fluxes on the interfaces
    nu_half=NaN(1,S.N);
    nu_half(2:S.N-1)=0.5*(nu(2:S.N-1)+nu(1:S.N-2));
    nu_half(1)=nu_left;
    nu_half(S.N)=nu_right;
    %Fluxes for nu equation
    nu_flux=zeros(S.N,1);
    nu_flux(2:S.N-1)=0.5*(u(1:S.N-2)+u(2:S.N-1));
    du_dh_left=-P_L(S.K,nu_left)*nu_left/zeta(nu_left);
    du_dh_right=-P_L(S.K,nu_right)*nu_right/zeta(nu_right);
    a=P_L(S.K,nu_right);
    b=zeta(nu_right);
    nu_flux(1)=-3*du_dh_left*S.h(1)/8+9*u(1)/8-u(2)/8;
    nu_flux(S.N)=3*du_dh_right*S.dh(S.N-1)/8-9*u(S.N-1)/8+u(S.N-2)/8;
    %Fluxes for u equation
    u_grad=2*(u(2:S.N-1)-u(1:S.N-2))./(S.dh(2:S.N-1)+S.dh(1:S.N-2));
    u_flux=zeros(1,S.N);
    u_flux(2:S.N-1)=a_1*P_L(S.K,nu_half(2:S.N-1))+(a_2*zeta(nu_half(2:S.N-1))./nu_half(2:S.N-1)).*u_grad;
    % The system of equations
    dnu_dt = nu_flux(2:S.N)-nu_flux(1:S.N-1);          
    du_dt =u_flux(2:S.N)'-u_flux(1:S.N-1)';     
    
    % Return the derivatives
    dydt = [dnu_dt; du_dt];
end

function dudh = du_dh(y,S)
dudh=NaN(S.N,1);
rho=y(1:S.N-1);
u=y(S.N:end);
dudh(2:S.N-1)=(u(2:S.N-1)-u(1:S.N-2))./(0.5*(S.dh(1:S.N-2)+S.dh(2:S.N-1)));
end
function y=P_L(K,nu)
y=K./nu;
end

function y=zeta(nu)
y=(2/3)*(nu).^(-2)./(nu-1);
end

function x=compute_grid(rho,rho_0,X)

%This function computes the new grid at each time step
%The equation for the new grid is given by: x_X=rho_0/rho, so we solve this
end

[Torsten]

gradients in u should give rise to a RHS for the first equation and therefore evolves nu, but it doesn't.

Where do you see in the code above that nu doesn't evolve ?

What is the original equation before you wrote it as a system of two equations ? What are the boundary conditions ?

[Mat]

The boundary conditions are: at both ends of the domain. The initial conditions are u(0,h)=0, and .

The two equations are the original set, but I integrate them to get the following set of ODEs:

I then approximate the derivatives and use the boundary conditions to find u on the boundary.

[Torsten]

When computing dnu_dt, you don't divide by S.dh in the inner grid points. I don't think this is correct.

[Mat]

You know, I think you might be right...

[Torsten]

It's hard to understand your discretization.

I'm also not sure whether you need a boundary value for u to be used in the first equation dnu/dt = du/dh.

It's a complicated system.

[Mat]

It is a complicated system, it's a mixed hyperbolic-parabolic system.

I need a boundary condition for u to get the end points for u, it's a Neumann to Diriclet map basically. I have included the discritisation.