% plot cubic function

u=[-0.1:0.01:1.1];
f=-u.*(u-a).*(u-1)+c;
plot(u,f)
hold on
axis([-0.1,1.1,-0.05,0.15]);
xlabel('u','fontsize',14);
ylabel('f(u)','fontsize',14);
subplot(211)
plot([-0.1,1.1],[0,0],'k:')  

% solve u'=f(u)

tspan=[0 10];  solve 0DE for t=0 to 10
u0=0.
[t1,u1]=ode45(@fcubic,tpsan,u0);