classdef fs2000_repl_solver2 < handle
	properties
		val
	end

	methods
		function this = fs2000_repl_solver2(k)
			this.val = this.runSolver(k);
		end

		function x = runSolver(this,k)
			fun = @this.steadyState;
			x0 = ones(1,17)*k
			opts=optimoptions('fsolve','Display','off','FunValCheck','off','MaxFunEvals',3400,'Algorithm','levenberg-marquardt');
			x = fsolve(fun,x0,opts);
		end
	end

	methods(Static)
		function F=steadyState(x)
			alp = 0.33;
			bet = 0.99;
			gam = 0.003;
			mst = 1.011;
			rho = 0.7;
			phi = 0.787;
			del = 0.02;

			da 			= x(1);
			gst 		= x(2);
			m 			= x(3);
			khst		= x(4);
			xist		= x(5);
			nust		= x(6);
			n			= x(7);
			p			= x(8);
			k			= x(9);
			l			= x(10);
			c			= x(11);
			d			= x(12);
			y			= x(13);
			r			= x(14);
			w			= x(15);
			gp_obs		= x(16);
			gy_obs		= x(17);
 
			F(1)		= -da + exp(gam); 
			F(2)		= -gst + 1/da;
			F(3)		= -m + mst;
			F(4)		= -khst + ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
			F(5)		= -xist + ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
			F(6)		= -nust + phi*mst^2/( (1-alp)*(1-phi)*bet*gst^alp*khst^alp );
			F(7)		= -n  + xist/(nust+xist);
			F(8)		= -p + xist + nust;
			F(9)		= -k + khst*n;
			F(10)		= -l + phi*mst*n/( (1-phi)*(1-n) );
			F(11)		= -c + mst/p;
			F(12)		= -d + l - mst + 1;
			F(13)		= -y + k^alp*n^(1-alp)*gst^alp;
			F(14)		= -r + mst/bet;
			F(15)		= -w + l/n;
			F(16)		= -gp_obs + m/da;
			F(17)		= -gy_obs + da;
		end
	end
end