Tolerance of ODE15s becomes too small

I want to solve the equations:

with boundary conditions , the initial conditions are .

I decided to solve this using the finite volume method, so that's why I don't need a boundary condition for the specific volume. I am told that this is a "stiff system", hence I am using ode15s. Do I need to fiddle with the ode15s options to get this thing to work?

%this code solves the following system of equations:
%v_t=u_h
%u_t=(zeta/v*u_h)_h
%dX/dt=u
%h is defined by h_x=1/v, and X is the new position of the interfacial
%points
%The boundary conditions for u is u(t,0)=0, u_h(t,1)=-v(t,1)
%The BC's for v, can be derived from ones on u. 
%here u is the velocity, and v is the specific volume defined as 1/rho, where rho is the density 
%The system is in conservative form, and the finite volume method is the best method for these types of equations 
clear;clc;
S=struct;
%---define physical constants 
N=200; %This is the number of cells in the simulation
S.N=N;
S.alpha=0;
c_p=1; %heat capacity.
rho_half(1:N)=0.8; %This is the initial density
eta_0=80; 
S.eta=eta_0;
S.nu_m_int=zeros(N+1,1);


%---Set up geometry
x=linspace(0,1,N+1); %Initial spacing of the cell interfaces and ends.
dx=x(2); %initial lengths of cells
L(1)=1; %Initial length

%---The solution of the system will be done via Newton-Raphson. The two
%variables v and u will be put into one vector y=[v u]' 
nu(1,:)=1./rho_half; %Initial specific volume
%This interpolates the density from the cell centres to cell boundaries.
rho_0(2:N)=0.5*(rho_half(1,N-1)+rho_half(2:N));
rho_0(1)=1.5*rho_half(1)-0.5*rho_half(2);
rho_0(N+1)=1.5*rho_half(N)-0.5*rho_half(N-1);
h=cumtrapz(x,rho_0); %calculates h, which is used throughout the timesteps.
nu_m(1:N)=1; %This is the matal powder density, the maximum possible density.
S.nu_m=nu_m;
dh=diff(h); %This is used to approximate the derivative
S.dh=dh;
S.h=h;
IC=zeros(2*N,1);
IC(1:N)=nu(1,:)';
%tspan = [0 5];
tspan = [0 1.8];
[t,y] = ode15s(@(t,y) sintering_RHS(y,S),tspan,IC);
figure(1)
plot(t,y(:,180))

figure(2)
plot(t,y(:,380))

function y=flux(u,nu,S)
%This computes the fluxes required for the finite volume method 
N=S.N;
dh=S.dh;
k=40;
a=1;
nu_face=0.5*(nu(1:N-1)+nu(2:N)); %This is specific volume on the cell interface
nu_L=1.5*nu(2)-0.5*nu(2);
nu_m=S.nu_m;
nu_m_int=S.nu_m_int;
nu_m_int(2:N)=0.5*(nu_m(1:N-1)+nu_m(2:N));
nu_m_int(N+1)=1.5*nu_m(N)-0.5*nu_m(N-1);
nu_m_int(1)=1.5*nu_m(1)-0.5*nu_m(2);
nu_end=1.5*nu(N)-0.5*nu(N-1);
y=zeros(2,N+1); %This vector will store the fluxes for the mass and momentum
y(1,2:N)=0.5*(u(2:N)+u(1:N-1)); %Fluxes for the mass
y(1,N+1)=u(N)-0.5*k*nu_end*dh(N);
y(2,2:N)=P_L(a,nu_face,nu_m_int(2:N))+(2*zeta(nu_face,nu_m_int(2:N),S)'./nu_face').*((u(2:N)'-u(1:N-1)')./(dh(2:N)+dh(1:N-1))); %The flux for the momentum equation
y(2,1)=P_L(a,nu_L,nu_m_int(1))+(zeta(nu_end,nu_m_int(1),S)/nu_end)*u(1); %the fluxes at the end
y(2,N+1)=-zeta(nu_end,nu_m_int(N+1),S);
end

function y=P_L(a,nu,nu_0) %This is the laplace pressure 
    x=1-nu_0./nu; %This is the porosity for an imcompressible metal powder
    y=a*(1-x').^2;
end

function y=zeta(z,nu_0,S)
eta_0=S.eta;
  x=1-nu_0./z;
  y = 2*(1-x).^3*eta_0./(3*x);
end

function f=sintering_RHS(y,S)
%the vector f has the stencil for each cell for both the mass and momentum
N=floor(length(y)/2);
f=zeros(1,2*N);
flux_1=flux(y(N+1:2*N),y(1:N),S); %This is flux at timestep i 
f_mass=flux_1(1,:);
f_mom=flux_1(2,:);
f(1:N)=f_mass(2:N+1)-f_mass(1:N); %These are the stencils coming from the system of equations 
f(N+1:2*N)=f_mom(2:N+1)-f_mom(1:N);
f=f';
end

[Torsten]

I decided to solve this using the finite volume method, so that's why I don't need a boundary condition for the specific volume.

The choice of the solution method cannot change the number of boundary conditions necessary to fix a solution for your problem.

Does your solution look reasonable up to the time when ode15s quits (i.e. up to 1.82 s) ?

[Mat Hunt]

No. I would have expected more diffusion, but I'm seeeing none.

[Torsten]

Did you consult some literature concerning the discretization for this type of PDE ? The system looks quite difficult.

[Mat Hunt]

I did. I also got tried to use my own Euler routine, and that worked very well but the time step was very small. I used a fully implicit scheme, and although I could get that working, it gives wonky results.

[Torsten]

With the initial and boundary conditions you wrote above, every numerical scheme should give you

u(h,t) = 0 for all h and all t
nu(h,t) = nu0(h) for all t    

as solution of your PDE system.

[Mat Hunt]

That's not what I get with my other numerical methods. The trivial solution is always an option.

[Torsten]

Starting from u(t=0,h) = 0 and u(t,h=0) = 0 together with du/dh(h=1) = 0, your scheme must produce u = 0 for all h and t. Otherwise, something is wrong with your discretization scheme or you implicitly use another boundary condition you are not aware of.

[Mat Hunt]

Thank you for your help. I'll have to sort this out myself it seems.

[Torsten]

Do you have a link to a publicly available article where the discretization of your PDE system is explained and where the equations have been solved ?

[Mat Hunt]

All I did was integrate over a cell to get the spacial discretisation.

[Torsten]

Special discretization schemes have been developed for special classes of equations (especially hyperbolic PDEs) in order to produce correct and stable solutions. Just using a standard approximation of the flux is often not good enough. That's why I asked if you consulted the literature for your system.

[Mat Hunt]

I took this method from LeVeque's book on finite volume methods.

[Torsten]

As far as I know, the book only treats scheme for hyperbolic systems of the form

du/dt + d(f(u))/dx = s(u,x,t)

I can't remember having seen schemes for 2nd order PDEs.

Your system is really special - a mixture of first and second order PDEs. I don't know either how to discretize it appropriately.