diff --git a/src/World3/population/common_pop/subsystems.jl b/src/World3/population/common_pop/subsystems.jl
index 3f38ee8..42bc2ea 100644
--- a/src/World3/population/common_pop/subsystems.jl
+++ b/src/World3/population/common_pop/subsystems.jl
@@ -10,7 +10,10 @@ function death_rate(; name, params=_params, inits=_inits, tables=_tables, ranges
     @variables dr(t) pop(t)
     @variables fpc(t) sopc(t) iopc(t) ppolx(t)
     @variables ehspc(t) = inits[:hsapc] # Line 22 Appendix A
-    @variables cdr(t) le(t) lmf(t) hsapc(t) lmhs(t) lmhs1(t) lmhs2(t) fpu(t) cmi(t) lmc(t) lmp(t)
+    @variables fpu(t) = inits[:fpu] 
+    @variables lmf(t) = inits[:lmf] 
+    @variables cmi(t) = inits[:cmi]
+    @variables cdr(t) le(t)  hsapc(t) lmhs(t) lmhs1(t) lmhs2(t)  lmc(t) lmp(t)
 
     eqs = [
         cdr ~ 1000.0 * dr / pop # Line 18 Appendix A
diff --git a/src/World3/population/pop15/initialisations.jl b/src/World3/population/pop15/initialisations.jl
index 8daac6f..e4ab326 100644
--- a/src/World3/population/pop15/initialisations.jl
+++ b/src/World3/population/pop15/initialisations.jl
@@ -2,7 +2,7 @@ _inits = Dict{Symbol, Float64}(
     :p1 => 5.3e7,
     :p2 => 2.1e8,
     :p3 => 2.2e8,
-    :p4 => 1.4e8,
+    :p4 => 1.8e8,
     :p5 => 1.6e8,
     :p6 => 1.5e8,
     :p7 => 1.3e8,
diff --git a/src/World3/population/pop15/plots.jl b/src/World3/population/pop15/plots.jl
index ce38825..8f6c5d3 100644
--- a/src/World3/population/pop15/plots.jl
+++ b/src/World3/population/pop15/plots.jl
@@ -1,6 +1,7 @@
+using DifferentialEquations
 function historicalrunsolution()
     isdefined(@__MODULE__, :_solution_historicalrun) && return _solution_historicalrun
-    global _solution_historicalrun = solve(historicalrun(), (1900, 2100))
+    global _solution_historicalrun = WorldDynamics.solve(historicalrun(), (1900, 2100); solver =  Rodas5())
     return _solution_historicalrun
 end
 
@@ -98,7 +99,7 @@ function fig_87(; kwargs...)
     parameters_2_87[:lt2] = 1900
 
     system = historicalrun(params=parameters_2_87)
-    solution = solve(system, (1900, 2100))
+    solution = WorldDynamics.solve(system, (1900, 2100); solver =  Rodas5())
 
     return plotvariables(solution, (t, 1900, 2100), _variables_a(); title="Fig. 2.87c", kwargs...)
 end
@@ -117,7 +118,7 @@ function fig_88(; kwargs...)
     parameters_2_88[:cfood] = 2500
 
     system = historicalrun(params=parameters_2_88)
-    solution = solve(system, (1900, 2100))
+    solution = WorldDynamics.solve(system, (1900, 2100); solver =  Rodas5())
 
     return plotvariables(solution, (t, 1900, 2100), _variables_a(); title="Fig. 2.88c", kwargs...)
 end
@@ -132,7 +133,7 @@ function fig_89(; kwargs...)
     parameters_2_89[:lt2] = 1900
 
     system = historicalrun(params=parameters_2_89)
-    solution = solve(system, (1900, 2100))
+    solution = WorldDynamics.solve(system, (1900, 2100); solver =  Rodas5())
 
     return plotvariables(solution, (t, 1900, 2100), _variables_a(); title="Fig. 2.89c", kwargs...)
 end
@@ -170,7 +171,7 @@ function fig94solution()
     new_equations[96] = pop.extra ~ 0.1 * pop.p6 + 0.25 * pop.p7 + 0.3 * pop.p8 + 0.25 * pop.p9 + 0.1 * pop.p10
 
     @named new_system = ODESystem(new_equations)
-    global _solution_94 = solve(new_system, (1900, 2100))
+    global _solution_94 = WorldDynamics.solve(new_system, (1900, 2100); solver =  Rodas5())
 
     return _solution_94
 end
@@ -207,7 +208,7 @@ function fig_96(; kwargs...)
     parameters_2_96[:fcest] = 1975
 
     system = historicalrun(params=parameters_2_96)
-    solution = solve(system, (1900, 2100))
+    solution = WorldDynamics.solve(system, (1900, 2100); solver =  Rodas5())
 
     return plotvariables(solution, (t, 1900, 2100), _variables_a(); title="Fig. 2.96c", kwargs...)
 end
@@ -223,7 +224,7 @@ function fig_97(; kwargs...)
     parameters_2_97[:zpgt] = 1975
 
     system = historicalrun(params=parameters_2_97)
-    solution = solve(system, (1900, 2100))
+    solution = WorldDynamics.solve(system, (1900, 2100); solver =  Rodas5())
 
     return plotvariables(solution, (t, 1900, 2100), _variables_a(); title="Fig. 2.97c", kwargs...)
 end
@@ -238,7 +239,7 @@ function fig_98(; kwargs...)
     parameters_2_98[:lt] = 2000
 
     system = historicalrun(params=parameters_2_98)
-    solution = solve(system, (1900, 2100))
+    solution = WorldDynamics.solve(system, (1900, 2100); solver =  Rodas5())
 
     return plotvariables(solution, (t, 1900, 2100), _variables_a(); title="Fig. 2.98c", kwargs...)
 end
@@ -254,7 +255,7 @@ function fig_99(; kwargs...)
     parameters_2_99[:fcest] = 1975
 
     system = historicalrun(params=parameters_2_99)
-    solution = solve(system, (1900, 2100))
+    solution = WorldDynamics.solve(system, (1900, 2100); solver =  Rodas5())
 
     return plotvariables(solution, (t, 1900, 2100), _variables_a(); title="Fig. 2.99c", kwargs...)
 end
@@ -271,7 +272,7 @@ function fig_100(; kwargs...)
     parameters_2_100[:zpgt] = 1975
 
     system = historicalrun(params=parameters_2_100)
-    solution = solve(system, (1900, 2100))
+    solution = WorldDynamics.solve(system, (1900, 2100); solver =  Rodas5())
 
     return plotvariables(solution, (t, 1900, 2100), _variables_a(); title="Fig. 2.100c", kwargs...)
 end