Why not just use some basic algebra?
I'll do it in MATLAB, as it will be faster that way.
By the way, using eps as a variable is just a bad idea, since eps is far too useful a tool on its own.
As1 = 510; fy=420; h = 650; d = 590; Es = 200000; fc = 21; As2 = 520; b = 350; beta1 = .85; T = As1*fy/1000;%kN
% I'll define c as symbolic here. syms c
E = (c-d)*.003./c; fs = Es*E; Cs = (fs-.85*fc)*As2/1000; Cc = beta1*fc*b*.85*c/1000; C = Cs + Cc; cond = C - T;
cond cond = (42483*c)/8000 + (104000*((3*c)/1000 - 177/100))/c - 111741/500
pretty(cond) / 3 c 177 \ | ---- - --- | 104000 42483 c \ 1000 100 / 111741 ------- + --------------------- - ------ 8000 c 500
Perhaps it is easier to just plot cond.
Note the value of h, since you require that c be between 0 and h.
h h = 650
So what are you doing with cond?
Pc = find(cond>=0); %STEP 3 Pc = min(Pc);
That is poor MATLAB. READ THE HELP FOR FIND, even though you should not be using use find here anyway.
If you want to find the first instance that makes your condition true, just do this:
Pc = find(cond>=0,1,'first');
Regardless, lets plot cond, as a function of c.
ezplot(cond[0,650]) grid on
The curve is a simple hyperbola. You want to find the smallest value of c that makes cond>=0.
Yes, you say in your question in one place that you want to find the LAST value of c. But that is inconsistent with your code, where you found the FIRST (SMALLEST) value of c.
So here is your function:
condfun = @(c) (42483*c)/8000 + (104000*((3*c)/1000 - 177/100))/c - 111741/500;
Now just use fzero on it. It crosses the axis at the point where cond==0.
[c0,fval] = fzero(condfun,[0.00001,650]) c0 = 178.04 fval = -1.4211e-13
Note that condfun is just slightly less than zero by an amount that is less than the tolerance for convergence. If that bothers you, then you can just bump c0 up by a tiny amount.
condfun(c0 + 1.e-12) ans = 1.1227e-11