%% Aufgabe 1a) und 1b)

%% Gleichgewichtszustände [a]

clear all;
close all;
clc;

syms v x1 k i x2 g m q p 

s = p*x1*x2;              % Symbolische Definition von s
b =(m*p*x1-q)/(m*p*k);    % Symbolische Definition von b

xp1 = v*x1*(1-x1/k)-s;    % Symbolische Definition von d/dx x1...3
xp2 = i*x2*(1-s/g)*b;
xp3 = m*s-q*x2;

f = [xp1;xp2;xp3];        % Alle drei Funktionen zu einem Vektor vereint
sol = solve(f==0,x1,x2);  % Ableitungen = 0 und umstellen nach x1, x2
a = [sol.x1,sol.x2];      % Lösungsteile für x1 und x2 zusammensetzen zu a

e = 0.05;     %Erosionslimit in Ressourceneinheiten [normiert]
g = 1;        %jährliches Produktionsziel [Ressourceneinheiten/Jahr]
k = 1.0;      %maximale Ressourcenkapazität [Ressourceneinheit]
m = 1.0;      %Preis für Ressource [Geldeinheit/Ressourceneinheit]
r = 0.1;      %Ressourcengenerationsrate [1/Jahr]

p = 1;        %Produktionsrate [1/(Kapital*Jahr)]
%p = 0.5;     %Produktionsrate [1/(Kapital*Jahr)]
q = 0.1;      %Kapital-Arbeitskosten [Geldeinheiten/(Kapital*Jahr)]
%q = 0.25;    %Kapital-Arbeitskosten [Geldeinheiten/(Kapital*Jahr)]
i = 0.1;      %Kapitalinvestitionsrate [1/Jahr]
%i = 0.25;    %Kapitalinvestitionsrate [1/Jahr]

n=1;


while n<=size(a,1)          % Substitution der symbolischen Var. mit konkreten Werten
    
    a(n,1) = subs(a(n,1));  % Lösungsteil x1 von a wird substituiert (und somit der konkrete Wert für x1 bestimmt)
    
    % Wert für x2 kann aufgrund Unkenntnis v noch nicht berechnet werden,
    % daher Bestimmung von v in if-Schleife:
    
    if a(n,1) < e,          % Welcher Wert wird v zugewiesen?
        v=0;                % Ist x1(t)<e wird v=0    
    else
        v=r;                % in allen anderen fällen (e>=0) wird v=r
    end
    
    a(n,2) = subs(a(n,2));  % Lösungsteil x2 von a wird substituiert
    
    n=n+1;
    
end

a