clc;
clear all;
warning off;
a = [.004; .006; .009];
b = [5.3; 5.5; 5.8]; 
c = [500; 400; 200];
PD=800;
DelP = 10;      % Error in Delta P is set to a high value
lambda = input('Enter estimated value of Lambda = ');
fprintf('\n ')
disp(['   Lambda      P1       P2        P3       DeltaP      DeltaLambda'])
iter = 0;      
while abs(DelP)  >= 0.001               
iter = iter + 1;                          
P = (lambda - b)./(2*a);
DelP =PD - sum(P);                        % Delta P
J = sum( ones(length(a), 1)./(2*a));      % sum(1/2*a)
Delambda = DelP/J;                        
disp([lambda, P(1), P(2), P(3), DelP, Delambda])
lambda = lambda + Delambda;             
end
lambda 
P
totalcost = sum(c + b.*P + a.*P.^2)
