Solving Linear Implicit Third Order Differential

Hey everyone,

I'm trying to model the viscoelastic behaviour of cells (with three relaxation times using Maxwell model) using the code below, which outputs an equation relating the indentation depth in an AFM test to the properties of the cells

clc
clear 
%housekeeping
syms h(t) A k v V G R y t a1(t) g1 g2 g3 tor1 tor2 tor3 h(B) B h0 h1 h2 h3 omega

%input symbols
L0 = (h(t)*((A*k)/(G)+ h(t)^0.5));
%LHS
L1 = diff((h(t)*((A*k)/(G)+ h(t)^0.5)),t);
L1 = simplify(L1);
%first deriviative LHS
L2 = diff((h(t)*((A*k)/(G)+ h(t)^0.5)),t,2);
L2 = simplify(L2);
%second deriviative LHS
L3 = diff((h(t)*((A*k)/(G)+ h(t)^0.5)),t,3);
L3 = simplify(L3);
%third deriviative LHS

H1 = int((h(B)*sym(exp(-(B-t)/tor1)))...
    ,B,0,t);
H2 = int((h(B)*sym(exp(-(B-t)/tor2)))...
    ,B,0,t);
H3 = int((h(B)*sym(exp(-(B-t)/tor3)))...
    ,B,0,t);
%set the integrals as constants

r0 = a1(t)-((A*k)/(g1*tor1)*H1)-((A*k)/(g2*tor2)*H2)-((A*k)/(g3*tor3)*H3);
%RHS(ramp)
r1 = diff(r0,t);
%first deriviative RHS(ramp)
r1 = simplify(r1); %simplify else it won't substitute properly later on
r2 = diff(r0,t,2);
r2 = simplify(r2);
%second deriviative RHS(ramp)
r3 = diff(r0,t,3);
r3 = collect(r3, tor3); %allows the substitution later
r3 = simplify(r3);
%third deriviative RHS(ramp)

%substitute the integrals out for a symbol
r0 = subs(r0, H1, h1);r0 = subs(r0, H2, h2);r0 = subs(r0, H3, h3);
r1 = subs(r1, H1, h1);r1 = subs(r1, H2, h2);r1 = subs(r1, H3, h3);
r2 = subs(r2, H1, h1);r2 = subs(r2, H2, h2);r2 = subs(r2, H3, h3);
r3 = subs(r3, H1, h1);r3 = subs(r3, H2, h2);r3 = subs(r3, H3, h3);

eqn0 = L0 == r0;
eqn1 = L1 == r1;
eqn2 = L2 == r2;
eqn3 = L3 == r3; %set up simultaneous equations between LHS and RHS

%begin substitution. Start at eqn3 and work backwards
H3 = solve(eqn3, h3); %H3 in terms of h1 and h2.
eqn2 = subs(eqn2, h3, H3); %removes h3 from circulation
eqn2 = simplify(eqn2);
H2 = solve(eqn2, h2); %H2 in terms of h1.
eqn1 = subs(eqn1, h3, H3);eqn1 = subs(eqn1, h2, H2);%removes h2 and h3 from circulation
eqn1 = simplify(eqn1);
H1 = solve(eqn1, h1);%H1 in terms of system variables
eqn0 = subs(eqn0, h3, H3);eqn0 = subs(eqn0, h2, H2);eqn0 = subs(eqn0, h1, H1);%removes h1, h2 and h3 from circulation 
%and gives the final equation
eqn0 = simplify(eqn0);

eqn0 = subs(eqn0, A, 3*(1-v)/(8*R^0.5)); %substitute A back in
eqn0 = subs(eqn0, a1(t),(A*k/G)*(1/G + 1/g1 + 1/g2 + 1/g3)*(V*t+h0)...
    - (A*k/G)*((V*tor1/g1)*(1-sym(exp(-t/tor1)) + (h0/g1)*exp(-t/tor1)))...
    - (A*k/G)*((V*tor2/g2)*(1-sym(exp(-t/tor2) + (h0/g2)*exp(-t/tor2))))...
    - (A*k/G)*((V*tor3/g3)*(1-sym(exp(-t/tor3) + (h0/g3)*exp(-t/tor3))))); 
%substitute a2(t) back in

Which outputs a third order implicit differential equation. I know I need to split it into three first order equations and solve it with an ode solver but I'm having a LOT of trouble. I've tried reduceDifferentialOrder/massMatrixForm/odeFunction and odeToVectorField/matlabFunction and daeFunction but I get a lot of error messages and I'm going a bit mad. This is what I've got right now:

[ramp,vars] = reduceDifferentialOrder(eqnr,h(t));
index = isLowIndexDAE(ramp,vars);
pramp = symvar(ramp);
pvars = symvar(vars);
params = setdiff(ramp,vars);
display(params)

f = daeFunction(ramp, vars,G, g1, g2, g3, h0, tor1, tor2, tor3);

G = 5E-2;
g1 = 5E-3; 
g2 = 5E-3;
g3 = 5E-3;
tor1 = 10;
tor2 = 100;
tor3 =1000;
y = 1; %ramp time
h0 = 0; 

F = @(t, Y, YP) f(t, Y, YP, G, g1, g2, g3, h0, tor1, tor2, tor3);

y0e = [0;V;0];
yp0e = zeros(3,1);
opt = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));
[y0, yp0] = decic(F, 0, y0e, [1;0;0], yp0e, [], opt);

[tF,yF] = ode15i(F, [0, y], y0, yp0, opt);

My main issue is that I get the error 'Index may be greater than one' when I run the decic code. If I skip that part, the integration fails. I don't know how to find the initial conditions manually. Please help!