Evaluating the Jacobian in a loop

Hi!

I have the following code providing me with 'steady state solutions' (pi*,c_sol*,z*):

    beta = 0.99 ; v = 21 ; gamma = 350 ; eps = 1 ; g = 0.2 ; sigma = 1; alpha = 0.7;
    c_sol = zeros(1,50);
    p = linspace(0.8,1.2,50);
    for   i = 1 : numel(p)
    pi_actual = p(i);
    fun = @(c)  (1 - beta)*pi_actual*(pi_actual-1) - (v/(alpha*gamma))*(c+g).^((1+eps)/alpha) - ((1-v)/gamma)*(c+g)*c.^(-sigma);
    c_sol(i) = fzero(fun,0.5);
    end
    plot(p,c_sol)
    z = p/beta;

Then I have the Jacobian as follows:

    syms c p R c2 p2
    outp = (c2*p2)/(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100);
    infl = very large equation
    J = [jacobian(outp, c2), jacobian(outp, p2) ; jacobian(infl, c2), jacobian(infl, p2)];

Now I want to calculate the eigenvalues for every point in the first script and plot the solutions versus pi*. Note:

    c = c_sol ; c2 = c_sol ; p2 = p ; R = z;

How can I make this a successful code? Thanks!