Uday, use one of the ode solvers:
function myODE()
IC = [0 0 1];
[t,X] = ode45(@ODEf,[0 10],IC);
f = X(:,1);
plot(t,f,'r')
xlabel('t')
ylabel('f')
grid
end
function df = ODEf(t, f)
df = zeros(3,1);
f1 = f(1);
f2 = f(2);
f3 = f(3);
df = [ f2;...
f3;...
-0.5*f1*f3];
end
