alternative for subs function

actually i am trying to solve matrix calculations with symbols when i use subs it does not evaluate but keeps it in huge number fractions huge means in millions and millions and i have to repeat for 10000 time with different values every time i get in previous step so please tell me some alternative i used simple but its not working and the whole process is so slow that in the 6th step it takes about 6 mins to solve my code is as follows:

I am attaching the m file also please help me out its urgent

thanks in advance

    clc
clear
format short
syms z1 z2 z3 z4 u z
z01 = input ('z01');
z02 = input ('z02');
z03 = input ('z03');
z04 = input ('z04');
u01 = input ('u01');
E = input ('E');
z0 = [z01;z02;z03;z04;u01;E];
z_dot = state()
J = sub(state())
s = 0.0001;
z_old = z0(2,1)
dJ = (der(z0(6,1)))
for i = 1:6
        z_new = (z_old) - s*(subs(dJ,{z1,z2,z3,z4,u},{z0(1,1),z_old,z0(3,1),z0(4,1),z0(5,1)}))
    z_old=z_new;
    hold on
    plot(i,z_old,'c*')

end

and the functions i used are as follows:

    function z_dot = state()
syms z1 z2 z3 z4 u
m = 0.23;
M = 1.0731;
I = 0.0079;
l = 0.3302;
b = 5.4;
g = 9.81;
%z = [z10;z20;z30;z40;u0];
%z = subs(z,[z1,z2,z3,z4])
F = u;
z1_dot = z3;
z2_dot = z4;
z3_dot = ((I+m*l^2)*F)/((M+m)*(I+m*l^2)-m^2*l^2*(cos(z2))^2) + (m*l*(I+m*l^2)*sin(z2)*(z4)^2)/((M+m)*(I+m*l^2)-m^2*l^2*(cos(z2))^2) - (b*(I+m*l^2)*z3)/((M+m)*(I+m*l^2)-m^2*l^2*(cos(z2))^2) + (m^2*l^2*sin(z2)*cos(z2))/((M+m)*(I+m*l^2)-m^2*l^2*(cos(z2))^2);
z4_dot = (-m*g*l*(M+m)*sin(z2))/((M+m)*(I+m*l^2)-m^2*l^2*(cos(z2))^2) - (m*l*cos(z2)*F)/((M+m)*(I+m*l^2)-m^2*l^2*(cos(z2))^2) - (m^2*l^2*sin(z2)*cos(z2)*(z4)^2)/((M+m)*(I+m*l^2)-m^2*l^2*(cos(z2))^2) - (m*l*b*cos(z2)*z3)/((M+m)*(I+m*l^2)-m^2*l^2*(cos(z2))^2);
z_dot = [z1_dot;z2_dot;z3_dot;z4_dot];
%z_dot = simple(subs(z_dot,[z1,z2,z3,z4,u],[z(1,1),z(2,1),z(3,1),z(4,1),z(5,1)]));
z_dot = simple(z_dot);
z_dot = simple(z_dot);
end

function derivativ = der(E)
format short
syms z1 z2 z3 z4 u
%E = input('E')
%z10 = input('z10')
%z20 = input('z20')
%z30 = input('z30')
%z40 = input('z40')
%u0 = input('u0')
%E = 0.01;
J2 = sub(state());
J2 = (subs(J2,{z1,z2,z3,z4,u},{z1,z2+E,z3,z4,u}));

J1 = sub(state()); J1 = (subs(J1,{z1,z2,z3,z4,u},{z1,z2,z3,z4,u}));

derivativ = ((J2-J1)/E);

end

    function substitute = sub(zk)
substitute = (norm(zk))^2;
end

[Neel]

my initial inputs are

z01 = 0

z02 = 0.2

z03 = 0

z04 = 0

u01 = 0

E = 0.01

Answer 1

[Walter Roberson]

You need to expect that sort of thing to happen when you do iteration of functions in rational fractions. It is needed to express the exact solution.

If you do not need the exact solution, have a look at vpa()

[Neel]

thank you so much this solves my whole problem