y0 = 1;
tspan = [0 1];

[t,y] = ode45('simple_rhs', tspan, y0);

ytrue = 3*exp(t) - t.^2 - 2*t - 2;
plot(t, ytrue, 'b-', t, y, 'rx')