diff --git a/examples/example01 - static_inverse_solve_MASTU.ipynb b/examples/example01 - static_inverse_solve_MASTU.ipynb
index afdf7dcb..ef2a2cf5 100644
--- a/examples/example01 - static_inverse_solve_MASTU.ipynb	
+++ b/examples/example01 - static_inverse_solve_MASTU.ipynb	
@@ -151,14 +151,18 @@
     "\n",
     "We can now instatiate a profile object that contains the chosen parameterisation of the toroidal plasma current density $J_p$ (i.e. on right hand side of the GS equation). We can then set the paramters for the chosen current density profiles. \n",
     "\n",
-    "A number of commonly used profile parameterisations exist in FreeGSNKE, including:\n",
-    "- `ConstrainPaxisIp`\n",
-    "- `ConstrainBetapIp`\n",
-    "- `Fiesta_Topeol`\n",
-    "- `Lao85`\n",
-    "- `TensionSpline`\n",
-    "\n",
-    "In this notebook, we will make use of the `ConstrainPaxisIp` (and `ConstrainBetapIp`) profiles (see [Jeon (2015)](https://link.springer.com/article/10.3938/jkps.67.843)). Others will be utilised in later notebooks. If there is a profile parameterisation you require that does not exist, please do create an issue. \n",
+    "The following table indicates which parameterisations are currently available in FreeGSNKE and for which types of equilibrium simulation they have been tested to work:\n",
+    "\n",
+    "| Name | Static | Evolutive (linear) | Evolutive (nonlinear) |\n",
+    "|-------|--------|------------------|---------------------|\n",
+    "| `ConstrainPaxisIp` | ✅ | ✅ | ✅ |\n",
+    "| `ConstrainBetapIp` | ✅ | ✅ | ✅ |\n",
+    "| `FiestaTopeol` | ✅ | ✅ | ✅ |\n",
+    "| `Lao85` | ✅ | ✅ | ✅ |\n",
+    "| `TensionSpline` | ✅ | ❌ | ❌ |\n",
+    "| `GeneralPprimeFFprime` | ✅ | ❌ | ❌ |\n",
+    "\n",
+    "In this notebook, we will make use of the `ConstrainPaxisIp` (and `ConstrainBetapIp`) profiles (see [Jeon (2015)](https://link.springer.com/article/10.3938/jkps.67.843)). Others will be utilised in later notebooks. If there is a profile parameterisation you require that does not exist, please do create an issue.\n",
     "\n",
     "Both `ConstrainPaxisIp` and `ConstrainBetapIp` are parameterised as follows:\n",
     " $$J_{p}(\\psi, R, Z) = \\lambda\\big[ \\beta_{0} \\frac{R}{R_{0}} \\left( 1-\\tilde{\\psi}^{\\alpha_m} \\right)^{\\alpha_n} + (1-\\beta_{0}) \\frac{R_0}{R} \\left( 1-\\tilde{\\psi}^{\\alpha_m} \\right)^{\\alpha_n} \\big] \\quad (R,Z) \\in \\Omega_p,$$\n",
diff --git a/examples/example02 - static_forward_solve_MASTU.ipynb b/examples/example02 - static_forward_solve_MASTU.ipynb
index a7a24ceb..856a1ddf 100644
--- a/examples/example02 - static_forward_solve_MASTU.ipynb	
+++ b/examples/example02 - static_forward_solve_MASTU.ipynb	
@@ -599,6 +599,117 @@
     "print(f\"Original beta_0 = {profiles_topeol.Beta0} vs. Fitted from Lao85 = {beta_0}.\")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Supposing an explicit parameterisation of the profiles is unknown and you have access to data describing the profiles. The `GeneralPprimeFFprime` profiles class enables you to use the raw data. \n",
+    "\n",
+    "As always the plasma current density is given by:\n",
+    "\n",
+    " $$J_{p}(\\psi, R, Z) = \\lambda\\big[ \\frac{R}{R_{0}} p'(\\tilde{\\psi}) + \\frac{R_0}{R} \\frac{1}{\\mu_0} F F'(\\tilde{\\psi}) \\big] \\quad (R,Z) \\in \\Omega_p, $$\n",
+    "\n",
+    "where the pressure, $p(\\tilde{\\psi})$, and toroidal magnetic field, $F(\\tilde{\\psi})$, profiles are now given by data arrays. \n",
+    "\n",
+    "By passing these data arrays either directly as `pprime_data` or `p_data` (and also `ffprime_data` or `f_data`), FreeGSNKE will interpolate the data and use them directly in the plasma current density above. \n",
+    "\n",
+    "The required parameters are:\n",
+    "- `Ip` (total plasma current).\n",
+    "- `fvac` ($rB_{tor}$,  vacuum toroidal field strength).\n",
+    "- `psi_n` ($\\tilde{\\psi}$, normalised poloidal flux at which data arrays corresponding, should be increasing).\n",
+    "- `pprime_data` or `p_data` (array of $p'(\\tilde{\\psi})$ or $p(\\tilde{\\psi})$ values at the normalised poloidal flux values in `psi_n`).\n",
+    "- `ffprime_data` or `f_data` (array of $FF'(\\tilde{\\psi})$ or $F(\\tilde{\\psi})$ values at the normalised poloidal flux values in `psi_n`).\n",
+    "- `Ip_logic` (if False, `Ip` is not used, if True, `Ip` is used to normalise $J_p$ and find $\\lambda$)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from freegsnke.jtor_update import GeneralPprimeFFprime\n",
+    "\n",
+    "psi_n = np.linspace(0,1,11)\n",
+    "\n",
+    "pprime_data = np.array([208061.72441505, 204110.95531309, 194359.61041113, 179811.09692964,\n",
+    "       161053.35395914, 138602.79928382, 113000.64896365,  84893.04354249,\n",
+    "        55164.04720078,  25298.53358472,      0.        ])\n",
+    "\n",
+    "ffprime_data = np.array([0.58961295, 0.57841712, 0.5507834 , 0.50955529, 0.45639891,\n",
+    "       0.3927777 , 0.32022539, 0.24057302, 0.1563259 , 0.07169191,\n",
+    "       0.        ])\n",
+    "\n",
+    "p_data = np.array([126175.75656716, 105510.77120406,  85544.30596718,  66798.45790761,\n",
+    "        49722.32952766,  34710.90959473,  22106.93373   ,  12194.53794856,\n",
+    "         5182.93833357,   1168.13190975,      0.        ])\n",
+    "\n",
+    "f_data = np.array([0.98240673, 0.92086948, 0.85722645, 0.79283784, 0.72925299,\n",
+    "       0.66837871, 0.61261316, 0.56490234, 0.52855955, 0.50657733,\n",
+    "       0.5       ])\n",
+    "\n",
+    "\n",
+    "# the following three setups are equivalent\n",
+    "profiles_general = GeneralPprimeFFprime(\n",
+    "    eq=eq,\n",
+    "    Ip=6e5,\n",
+    "    fvac=0.5,\n",
+    "    psi_n=psi_n,\n",
+    "    pprime_data=pprime_data,\n",
+    "    ffprime_data=ffprime_data,\n",
+    "    p_data=p_data,\n",
+    "    f_data=f_data,\n",
+    ")\n",
+    "\n",
+    "# profiles_general = GeneralPprimeFFprime(\n",
+    "#     eq=eq,\n",
+    "#     Ip=6e5,\n",
+    "#     fvac=0.5,\n",
+    "#     psi_n=psi_n,\n",
+    "#     pprime_data=pprime_data,\n",
+    "#     ffprime_data=ffprime_data,\n",
+    "#     p_data=None,\n",
+    "#     f_data=None,\n",
+    "# )\n",
+    "\n",
+    "# profiles_general = GeneralPprimeFFprime(\n",
+    "#     eq=eq,\n",
+    "#     Ip=6e5,\n",
+    "#     fvac=0.5,\n",
+    "#     psi_n=psi_n,\n",
+    "#     pprime_data=None,\n",
+    "#     ffprime_data=None,\n",
+    "#     p_data=p_data,\n",
+    "#     f_data=f_data,\n",
+    "# )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# instatiate new equilibrium object\n",
+    "eq_general = deepcopy(eq)\n",
+    "\n",
+    "# call solver with new profile object\n",
+    "GSStaticSolver.solve(eq=eq_general, \n",
+    "                     profiles=profiles_general, \n",
+    "                     constrain=None, \n",
+    "                     target_relative_tolerance=1e-9)\n",
+    "\n",
+    "\n",
+    "# plot the resulting equilbria \n",
+    "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n",
+    "ax1.grid(True, which='both')\n",
+    "eq_general.plot(axis=ax1, show=False)\n",
+    "eq_general.tokamak.plot(axis=ax1, show=False)\n",
+    "ax1.set_xlim(0.1, 2.15)\n",
+    "ax1.set_ylim(-2.25, 2.25)\n",
+    "plt.tight_layout()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
diff --git a/examples/example03 - extracting_equilibrium_quantites.ipynb b/examples/example03 - extracting_equilibrium_quantites.ipynb
index 5f54f242..7cfa0f9c 100644
--- a/examples/example03 - extracting_equilibrium_quantites.ipynb	
+++ b/examples/example03 - extracting_equilibrium_quantites.ipynb	
@@ -906,7 +906,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# plot the input p' and FF' profiles\n",
+    "# plot the p' and FF' profiles\n",
     "\n",
     "psi_n = eq.psiN_1D(N=65)\n",
     "\n",
@@ -929,15 +929,21 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# plot q profile\n",
+    "# plot the p and F profiles\n",
     "\n",
-    "psi_n = np.linspace(0.01,0.99,65)  # values of q at 0 and 1 can be problematic\n",
+    "psi_n = eq.psiN_1D(N=65)\n",
     "\n",
-    "fig1, ax1 = plt.subplots(1, 1, figsize=(6,6), dpi=80)\n",
+    "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,6), dpi=80)\n",
     "ax1.grid(zorder=0, alpha=0.75)\n",
-    "ax1.plot(psi_n, eq.q(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
+    "ax1.plot(psi_n, profiles.pressure(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
     "ax1.set_xlabel(r'$\\hat{\\psi}$')\n",
-    "ax1.set_ylabel(r\"$q(\\hat{\\psi})$\")\n"
+    "ax1.set_ylabel(r\"$p(\\hat{\\psi})$\")\n",
+    "ax1.ticklabel_format(axis='y', scilimits=(0,0))\n",
+    "\n",
+    "ax2.grid(zorder=0, alpha=0.75)\n",
+    "ax2.plot(psi_n, profiles.fpol(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
+    "ax2.set_xlabel(r'$\\hat{\\psi}$')\n",
+    "ax2.set_ylabel(r\"$F(\\hat{\\psi})$\")\n"
    ]
   },
   {
@@ -946,14 +952,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# plot fpol\n",
+    "# plot q profile\n",
+    "\n",
+    "psi_n = np.linspace(0.01,0.99,65)  # values of q at 0 and 1 can be problematic\n",
     "\n",
     "fig1, ax1 = plt.subplots(1, 1, figsize=(6,6), dpi=80)\n",
     "ax1.grid(zorder=0, alpha=0.75)\n",
-    "ax1.plot(psi_n, eq.fpol(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
+    "ax1.plot(psi_n, eq.q(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
     "ax1.set_xlabel(r'$\\hat{\\psi}$')\n",
-    "ax1.set_ylabel(r\"$fpol(\\hat{\\psi})$\")\n",
-    "\n"
+    "ax1.set_ylabel(r\"$q(\\hat{\\psi})$\")\n"
    ]
   },
   {
diff --git a/examples/example05 - evolutive_forward_solve.ipynb b/examples/example05 - evolutive_forward_solve.ipynb
index e9c3435e..d04a0efd 100644
--- a/examples/example05 - evolutive_forward_solve.ipynb	
+++ b/examples/example05 - evolutive_forward_solve.ipynb	
@@ -11,13 +11,11 @@
     "To do this, we need to:\n",
     "- build the tokamak machine.\n",
     "- instatiate a GS equilibrium (to be used as an initial condition for the evolutive solver).\n",
-    "- calculate a vertical instability growth rate for this equilibrium and carry out passive structure mode removal via a normal mode decomposition (i.e. removing modes that have little effect on the evolution). \n",
-    "- define time-dependent plasma current density profile parameters and coil voltages.\n",
+    "- calculate a vertical instability growth rate for this equilibrium and carry out passive structure mode removal via a normal mode decomposition (i.e. removing modes that have little effect on the evolution). More details on this in a later notebook. \n",
+    "- define time-dependent coil voltages, plasma current density profile parameters, and plasma resistivity.\n",
     "- evolve the active coil currents, the total plasma current, and the equilbirium using these profile parameters and voltages by solving the circuit equations alongside the GS equation.\n",
     "\n",
-    "Refer to the paper by [Amorisco et al. (2024)](https://pubs.aip.org/aip/pop/article/31/4/042517/3286904/FreeGSNKE-A-Python-based-dynamic-free-boundary) for more details. \n",
-    "\n",
-    "We should note that here we will use **fixed** (time-independent) profile parameters and voltages to simulate a VDE, however, truly time-dependent parameters would be required to simulate a plasma shot (see future notebooks). \n"
+    "Refer to the paper by [Amorisco et al. (2024)](https://pubs.aip.org/aip/pop/article/31/4/042517/3286904/FreeGSNKE-A-Python-based-dynamic-free-boundary) for more details."
    ]
   },
   {
@@ -37,7 +35,8 @@
     "import matplotlib.pyplot as plt\n",
     "from copy import deepcopy\n",
     "from IPython.display import display, clear_output\n",
-    "import pickle"
+    "import pickle\n",
+    "import copy"
    ]
   },
   {
@@ -220,11 +219,11 @@
     "- `min_dIy_dI`: threshold value below which passive structure normal modes are dropped. Modes with norm(d(Iy)/dI)<`min_dIy_dI` are dropped, which filters out modes that do not actually couple with the plasma.\n",
     "- `max_mode_frequency`: threshold value for characteristic frequencies above which passive structure normal modes are dropped (i.e. the fast modes).\n",
     "\n",
-    "Other customisable inputs are available, do see the documentation for more details. For example, one may explicitly set your own resistance and inductance matries for the tokamak machine. \n",
+    "Other customisable inputs are available, do see the documentation or the later notebook on \"Growth Rates\" for more details. For example, one may explicitly set your own resistance and inductance matries for the tokamak, rather than the geometrical value calculated internally in FreeGSNKE.\n",
     "\n",
     "The solver can be used on different equilibria and/or profiles, but these need to have the same machine, domain, and limiter as the one used at the time of the solver instantiation. For different machines, a new time-evolutive solver should be created.\n",
     "\n",
-    "The input equilibrium and profile functions are also used as the expansion point around which the dynamics are linearised."
+    "The input equilibrium and profile functions are also used as the expansion point around which the dynamics are linearised (when using the linear solver or calculating linear growth rates). "
    ]
   },
   {
@@ -245,44 +244,6 @@
     ")"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "#### Set active coil voltages\n",
-    "\n",
-    "In this example, we will evolve a plasma in absence of any control policy or current drive.\n",
-    "\n",
-    "Just as an example, the following calculates active voltages to be applied to the poloidal field coils (and Solenoid) using $V = RI$, with current values as defined by the initial equilibrium (i.e. we have constant voltages).\n",
-    "\n",
-    "In most FreeGSNKE use cases, these active voltages will be determined by a control policy."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "voltages = (stepping.vessel_currents_vec*stepping.evol_metal_curr.R)[:stepping.evol_metal_curr.n_active_coils] "
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "To start, the solver is prepared by setting the initial conditions."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "stepping.initialize_from_ICs(eq, profiles)"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -312,8 +273,27 @@
     "history_triangularity = [stepping.eq1.triangularity()]\n",
     "history_squareness = [stepping.eq1.squareness()[1]]\n",
     "history_area = [stepping.eq1.separatrix_area()]\n",
-    "history_length = [stepping.eq1.separatrix_length()]\n",
-    "\n"
+    "history_length = [stepping.eq1.separatrix_length()]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Set time-dependent active coil voltages, profile parameters, and plasma resistivity\n",
+    "\n",
+    "To simulate the VDE in the example, we will evolve a plasma in absence of any control policy or current drive. The following calculates active voltages to be applied to the poloidal field coils (and Solenoid) using $V = RI$, with current values as defined by the initial equilibrium (i.e. we have constant voltages over time). For more realistic plasma simulation, these active voltages would be truly time-dependent and determined by a control policy.\n",
+    "\n",
+    "In this example, the profile parameters and the plasma resistivity will also remain constant over time. In the next example, we will show how to vary these."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "voltages = (stepping.vessel_currents_vec*stepping.evol_metal_curr.R)[:stepping.evol_metal_curr.n_active_coils] "
    ]
   },
   {
@@ -323,7 +303,7 @@
     "#### Call the solver (linear)\n",
     "Finally, we call the time-evolutive solver with `stepping.nlstepper()` sequentially until we reach the preset end time.\n",
     "\n",
-    "The following demonstrates a solely linear evolution of the plasma by setting `linear_only=True`."
+    "The following demonstrates a solely linear evolution of the plasma by setting `linear_only=True`. This makes use of the Jacobians calculated when the nonlinear solver object was initialised. Note that the linear solution will only hold for equilibrium \"close to\" that which was used to initialise the nonlinear solver object. "
    ]
   },
   {
@@ -332,15 +312,20 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# initialise the solver with the initial equilibrium/profiles\n",
+    "stepping.initialize_from_ICs(eq, profiles)\n",
+    "\n",
     "# loop over time steps\n",
     "while counter<max_count:\n",
     "    print(f'Step: {counter}/{max_count-1}')\n",
     "    print(f'--- t = {t:.2e}')\n",
     "\n",
-    "    # carry out the time step\n",
+    "    # carry out the time step (profile parameters and resistivity are constant at each time step)\n",
     "    stepping.nlstepper(\n",
-    "        active_voltage_vec=voltages, \n",
+    "        active_voltage_vec=voltages,   # same voltages used at each time step\n",
     "        linear_only=True,\n",
+    "        verbose=False,\n",
+    "        max_solving_iterations=50,\n",
     "    )\n",
     "\n",
     "    # store information on the time step\n",
@@ -369,23 +354,6 @@
     "history_length = np.array(history_length)\n"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n",
-    "\n",
-    "ax1.grid(True, which='both')\n",
-    "history_equilibria[-1].plot(axis=ax1, show=False)\n",
-    "eq.tokamak.plot(axis=ax1, show=False)\n",
-    "ax1.set_xlim(0.1, 2.15)\n",
-    "ax1.set_ylim(-2.25, 2.25)\n",
-    "\n",
-    "plt.tight_layout()"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -585,9 +553,7 @@
    "source": [
     "## Limited plasma example\n",
     "\n",
-    "Finally, we examine a limiter configuration (there is no fundamental difference from the process described above).\n",
-    "\n",
-    "First we re-instantiate the equilibrium object."
+    "In this example, we examine a limiter configuration plasma and show how to use time-varying plasma profile parameters and plasma resisitivity."
    ]
   },
   {
@@ -596,40 +562,73 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "from freegsnke.jtor_update import ConstrainPaxisIp\n",
+    "from freegsnke import equilibrium_update, GSstaticsolver\n",
+    "\n",
+    "# re-instantiate the equilibrium object\n",
     "eq = equilibrium_update.Equilibrium(\n",
     "    tokamak=tokamak,\n",
     "    Rmin=0.1, Rmax=2.0,\n",
     "    Zmin=-2.2, Zmax=2.2,\n",
     "    nx=65, \n",
     "    ny=129\n",
-    ")  "
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We load a set of coil currents that will give us a limiter plasma in our current tokamak."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
+    ")  \n",
+    "\n",
+    "# load coil currents\n",
     "with open('data/simple_limited_currents_PaxisIp.pk', 'rb') as f:\n",
     "    current_values = pickle.load(f)\n",
     "\n",
     "for key in current_values.keys():\n",
-    "    eq.tokamak.set_coil_current(key, current_values[key])"
+    "    eq.tokamak.set_coil_current(key, current_values[key])\n",
+    "\n",
+    "\n",
+    "# define profiles\n",
+    "profiles = ConstrainPaxisIp(\n",
+    "    eq=eq,\n",
+    "    paxis=6e3,\n",
+    "    Ip=4e5,\n",
+    "    fvac=0.5,\n",
+    "    alpha_m=1.8,\n",
+    "    alpha_n=1.2\n",
+    ")\n",
+    "\n",
+    "# solve\n",
+    "GSStaticSolver = GSstaticsolver.NKGSsolver(eq)    \n",
+    "GSStaticSolver.solve(\n",
+    "    eq, \n",
+    "    profiles, \n",
+    "    target_relative_tolerance=1e-8, \n",
+    "    verbose=False\n",
+    "    )\n",
+    "\n",
+    "# visualise:\n",
+    "# the red dashed line is the flux surface through the first X-point, however, the actual last closed flux surface \n",
+    "# in limiter equilibria is displayed as a full black line. \n",
+    "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n",
+    "ax1.grid(True, which='both')\n",
+    "eq.plot(axis=ax1, show=False)\n",
+    "eq.tokamak.plot(axis=ax1, show=False)\n",
+    "ax1.set_xlim(0.1, 2.15)\n",
+    "ax1.set_ylim(-2.25, 2.25)\n",
+    "plt.tight_layout()"
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
+   "outputs": [],
    "source": [
-    "Now we set the profiles."
+    "from freegsnke import nonlinear_solve\n",
+    "\n",
+    "# initialise the nonlinear solver object\n",
+    "stepping = nonlinear_solve.nl_solver(\n",
+    "    eq=eq, \n",
+    "    profiles=profiles, \n",
+    "    full_timestep=.5e-3, \n",
+    "    plasma_resistivity=1e-6, \n",
+    "    GSStaticSolver=GSStaticSolver   \n",
+    ")"
    ]
   },
   {
@@ -638,22 +637,18 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "profiles = ConstrainPaxisIp(\n",
-    "    eq=eq,\n",
-    "    paxis=6e3,\n",
-    "    Ip=4e5,\n",
-    "    fvac=0.5,\n",
-    "    alpha_m=1.8,\n",
-    "    alpha_n=1.2\n",
-    ")"
+    "# number of time steps to simulate\n",
+    "max_count = 50"
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
+   "outputs": [],
    "source": [
-    "The final part of the setup is to do a forward solve and visualise the limiter equilibrium.\n",
-    "The red dashed line is the flux surface through the first X-point, however, the actual last closed flux surface in limiter equilibria is displayed as a full black line. "
+    "# determine the (constant) active voltages to be applied at each time step\n",
+    "voltages = (stepping.vessel_currents_vec*stepping.evol_metal_curr.R)[:stepping.evol_metal_curr.n_active_coils] "
    ]
   },
   {
@@ -662,22 +657,56 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "GSStaticSolver.solve(eq, profiles, target_relative_tolerance=1e-8, verbose=True)\n",
+    "# we define some time-dependent plasma current density profile parameters here\n",
+    "alpha_m = np.tile(profiles.alpha_m, max_count+1)\n",
+    "alpha_m -= (0.1 * np.sin(0.05 * np.pi * np.arange(max_count+1))) # we add some perturbation\n",
     "\n",
-    "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n",
-    "ax1.grid(True, which='both')\n",
-    "eq.plot(axis=ax1, show=False)\n",
-    "eq.tokamak.plot(axis=ax1, show=False)\n",
-    "ax1.set_xlim(0.1, 2.15)\n",
-    "ax1.set_ylim(-2.25, 2.25)\n",
-    "plt.tight_layout()"
+    "alpha_n = np.tile(profiles.alpha_n, max_count+1)\n",
+    "alpha_n += (0.1 * np.sin(0.1 * np.pi * np.arange(max_count+1))) # we add some perturbation\n",
+    "\n",
+    "paxis = np.tile(profiles.paxis, max_count+1)\n",
+    "paxis += (0.1 * np.sin(0.01 * np.pi * np.arange(max_count+1))) # we add some perturbation\n",
+    "\n",
+    "\n",
+    "# plot of profiles\n",
+    "fig, axs = plt.subplots(3, 1, figsize=(8, 6), dpi=80, constrained_layout=True)\n",
+    "axs_flat = axs.flat\n",
+    "\n",
+    "axs_flat[0].plot(alpha_m, label='alpha_m')\n",
+    "axs_flat[0].set_xticklabels([])\n",
+    "axs_flat[0].set_ylabel('Parameter value')\n",
+    "axs_flat[0].legend()\n",
+    "axs_flat[0].grid()\n",
+    "\n",
+    "axs_flat[1].plot(alpha_n, label='alpha_n')\n",
+    "axs_flat[1].set_xticklabels([])\n",
+    "axs_flat[1].set_ylabel('Parameter value')\n",
+    "axs_flat[1].legend()\n",
+    "axs_flat[1].grid()\n",
+    "\n",
+    "axs_flat[2].plot(paxis, label='paxis')\n",
+    "axs_flat[2].set_xlabel('Time Step')\n",
+    "axs_flat[2].set_ylabel('Parameter value')\n",
+    "axs_flat[2].legend()\n",
+    "axs_flat[2].grid()\n"
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
+   "outputs": [],
    "source": [
-    "Finally, we can do the nonlinear evolutive solve as before. First, we re-initialize the solver, which calculates the linearization of the dynamics relevant to this new equilibrium."
+    "# we also define a time-dependent plasma resisitivity. \n",
+    "resisitivity = np.tile(stepping.plasma_resistivity, max_count+1)\n",
+    "resisitivity += (1e-7 * np.sin(0.05 * np.pi * np.arange(max_count+1))) # we add some perturbation\n",
+    "\n",
+    "# plot of resisitivity evolution\n",
+    "fig, axs = plt.subplots(1, 1, figsize=(8, 3), dpi=80)\n",
+    "axs.plot(resisitivity)\n",
+    "axs.set_xlabel('Time Step')\n",
+    "axs.set_ylabel('Resisitivity')\n",
+    "axs.grid()\n"
    ]
   },
   {
@@ -686,14 +715,65 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from freegsnke import nonlinear_solve\n",
+    "# initialising some variables for iteration and logging\n",
+    "counter = 0\n",
+    "t = 0\n",
     "\n",
-    "stepping = nonlinear_solve.nl_solver(\n",
-    "    eq=eq, profiles=profiles, \n",
-    "    full_timestep=.5e-3, \n",
-    "    plasma_resistivity=1e-6, \n",
-    "    GSStaticSolver=GSStaticSolver   \n",
-    ")"
+    "# reset the solver object by resetting the initial conditions\n",
+    "stepping.initialize_from_ICs(eq, profiles)\n",
+    "\n",
+    "history_times = [t]\n",
+    "history_currents = [stepping.currents_vec]\n",
+    "history_equilibria = [deepcopy(stepping.eq1)]\n",
+    "history_o_points = [stepping.eq1.opt[0]]\n",
+    "history_elongation = [stepping.eq1.geometricElongation()]\n",
+    "history_triangularity = [stepping.eq1.triangularity()]\n",
+    "history_squareness = [stepping.eq1.squareness()[1]]\n",
+    "history_area = [stepping.eq1.separatrix_area()]\n",
+    "history_length = [stepping.eq1.separatrix_length()]\n",
+    "\n",
+    "# simulate using the linear solver\n",
+    "while counter<max_count:\n",
+    "\n",
+    "    print(f'Step: {counter}/{max_count-1}')\n",
+    "    print(f'--- t = {t:.2e}')\n",
+    "    \n",
+    "    # carry out the time step\n",
+    "    stepping.nlstepper(\n",
+    "        active_voltage_vec = voltages, # same voltages applied at each step\n",
+    "        linear_only = True,\n",
+    "        profiles_parameters = { # different at each step\n",
+    "            \"alpha_m\": alpha_m[counter],\n",
+    "            \"alpha_n\": alpha_n[counter],\n",
+    "            \"paxis\": paxis[counter],\n",
+    "            },\n",
+    "        plasma_resistivity = resisitivity[counter] # different at each step\n",
+    "        )\n",
+    "\n",
+    "    # store information on the time step\n",
+    "    t += stepping.dt_step\n",
+    "    history_times.append(t)\n",
+    "    counter += 1\n",
+    "    \n",
+    "    # store time-advanced equilibrium, currents, and profiles (+ other quantites of interest)\n",
+    "    history_currents.append(stepping.currents_vec)\n",
+    "    history_equilibria.append(deepcopy(stepping.eq1))\n",
+    "    history_o_points.append(stepping.eq1.opt[0])\n",
+    "    history_elongation.append(stepping.eq1.geometricElongation())\n",
+    "    history_triangularity.append(stepping.eq1.triangularity())\n",
+    "    history_squareness.append(stepping.eq1.squareness()[1])\n",
+    "    history_area.append(stepping.eq1.separatrix_area())\n",
+    "    history_length.append(stepping.eq1.separatrix_length())\n",
+    "\n",
+    "# transform lists to arrays\n",
+    "history_currents = np.array(history_currents)\n",
+    "history_times = np.array(history_times)\n",
+    "history_o_points = np.array(history_o_points)\n",
+    "history_elongation = np.array(history_elongation)\n",
+    "history_triangularity = np.array(history_triangularity)\n",
+    "history_squareness = np.array(history_squareness)\n",
+    "history_area = np.array(history_area)\n",
+    "history_length = np.array(history_length)"
    ]
   },
   {
@@ -702,19 +782,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# recalculate the active voltages using the new currents\n",
-    "U_active = (stepping.vessel_currents_vec*stepping.evol_metal_curr.R)[:stepping.evol_metal_curr.n_active_coils] \n",
-    "\n",
-    "# number of time steps to simulate\n",
-    "max_count = 50\n",
-    "\n",
-    "# reset the solver object by resetting the initial conditions\n",
-    "stepping.initialize_from_ICs(eq, profiles)\n",
-    "\n",
     "# initialising some variables for iteration and logging\n",
     "counter = 0\n",
     "t = 0\n",
     "\n",
+    "# reset the solver object by resetting the initial conditions\n",
+    "stepping.initialize_from_ICs(eq, profiles)\n",
+    "\n",
     "history_times_nl = [t]\n",
     "history_currents_nl = [stepping.currents_vec]\n",
     "history_equilibria_nl = [deepcopy(stepping.eq1)]\n",
@@ -725,7 +799,7 @@
     "history_area_nl = [stepping.eq1.separatrix_area()]\n",
     "history_length_nl = [stepping.eq1.separatrix_length()]\n",
     "\n",
-    "# Simulate the forward evolution\n",
+    "# simulate using the nonlinear solver\n",
     "while counter<max_count:\n",
     "\n",
     "    print(f'Step: {counter}/{max_count-1}')\n",
@@ -733,9 +807,15 @@
     "    \n",
     "    # carry out the time step\n",
     "    stepping.nlstepper(\n",
-    "        active_voltage_vec=U_active, \n",
-    "        linear_only=False,\n",
-    "    )\n",
+    "        active_voltage_vec = voltages, # same voltages applied at each step\n",
+    "        linear_only = False,\n",
+    "        profiles_parameters = { # different profile parameters at each step\n",
+    "            \"alpha_m\": alpha_m[counter],\n",
+    "            \"alpha_n\": alpha_n[counter],\n",
+    "            \"paxis\": paxis[counter],\n",
+    "            },\n",
+    "        plasma_resistivity = resisitivity[counter] # different resistivity at each step\n",
+    "        )\n",
     "\n",
     "    # store information on the time step\n",
     "    t += stepping.dt_step\n",
@@ -773,48 +853,57 @@
     "fig, axs = plt.subplots(3, 3, figsize=(15, 10), dpi=80, constrained_layout=True)\n",
     "axs_flat = axs.flat\n",
     "\n",
+    "axs_flat[0].plot(history_times, history_o_points[:, 0],'k+', label='linear')\n",
     "axs_flat[0].plot(history_times_nl, history_o_points_nl[:, 0],'rx', label='nonlinear')\n",
     "axs_flat[0].set_xlabel('Time')\n",
     "axs_flat[0].set_ylabel('Magnetic axis, $R$ [$m$]')\n",
     "axs_flat[0].legend()\n",
     "axs_flat[0].grid()\n",
     "\n",
+    "axs_flat[1].plot(history_times, abs(history_o_points[:, 1]),'k+')\n",
     "axs_flat[1].plot(history_times_nl, abs(history_o_points_nl[:, 1]),'rx')\n",
     "axs_flat[1].set_yscale('log')\n",
     "axs_flat[1].set_xlabel('Time')\n",
     "axs_flat[1].set_ylabel('abs(Magnetic axis, $Z$) [$m$]')\n",
     "axs_flat[1].grid()\n",
     "\n",
+    "axs_flat[2].plot(history_times, history_o_points[:, 2],'k+')\n",
     "axs_flat[2].plot(history_times_nl, history_o_points_nl[:, 2],'rx')\n",
     "axs_flat[2].set_xlabel('Time')\n",
     "axs_flat[2].set_ylabel('Magnetic axis, $\\psi$')\n",
     "axs_flat[2].grid()\n",
     "\n",
+    "axs_flat[3].plot(history_times, history_currents[:,-1]*stepping.plasma_norm_factor,'k+')\n",
     "axs_flat[3].plot(history_times_nl, history_currents_nl[:,-1]*stepping.plasma_norm_factor,'rx')\n",
     "axs_flat[3].set_xlabel('Time')\n",
     "axs_flat[3].set_ylabel('Plasma current')\n",
     "axs_flat[3].grid()\n",
     "\n",
+    "axs_flat[4].plot(history_times, history_elongation,'k+')\n",
     "axs_flat[4].plot(history_times_nl, history_elongation_nl,'rx')\n",
     "axs_flat[4].set_xlabel('Time')\n",
     "axs_flat[4].set_ylabel('Geometric elongation')\n",
     "axs_flat[4].grid()\n",
     "\n",
+    "axs_flat[5].plot(history_times, history_triangularity,'k+')\n",
     "axs_flat[5].plot(history_times_nl, history_triangularity_nl,'rx')\n",
     "axs_flat[5].set_xlabel('Time')\n",
     "axs_flat[5].set_ylabel('Triangularity')\n",
     "axs_flat[5].grid()\n",
     "\n",
+    "axs_flat[6].plot(history_times, history_squareness,'k+')\n",
     "axs_flat[6].plot(history_times_nl, history_squareness_nl,'rx')\n",
     "axs_flat[6].set_xlabel('Time')\n",
     "axs_flat[6].set_ylabel('Squarenss')\n",
     "axs_flat[6].grid()\n",
     "\n",
+    "axs_flat[7].plot(history_times, history_area,'k+')\n",
     "axs_flat[7].plot(history_times_nl, history_area_nl,'rx')\n",
     "axs_flat[7].set_xlabel('Time')\n",
     "axs_flat[7].set_ylabel('LCFS area')\n",
     "axs_flat[7].grid()\n",
     "\n",
+    "axs_flat[8].plot(history_times, history_length,'k+')\n",
     "axs_flat[8].plot(history_times_nl, history_length_nl,'rx')\n",
     "axs_flat[8].set_xlabel('Time')\n",
     "axs_flat[8].set_ylabel('LCFS circumference')\n",
@@ -831,21 +920,29 @@
     "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 8), dpi=80)\n",
     "\n",
     "ax1.grid(True, which='both')\n",
-    "history_equilibria_nl[-1].plot(axis=ax1, show=False)\n",
+    "history_equilibria[-1].plot(axis=ax1, show=False)\n",
     "eq.tokamak.plot(axis=ax1, show=False)\n",
     "ax1.set_xlim(0.1, 2.15)\n",
     "ax1.set_ylim(-2.25, 2.25)\n",
-    "ax1.set_title('t=0')\n",
+    "ax1.set_title(\"Linear solve\")\n",
     "\n",
     "ax2.grid(True, which='both')\n",
     "history_equilibria_nl[-1].plot(axis=ax2, show=False)\n",
     "eq.tokamak.plot(axis=ax2, show=False)\n",
     "ax2.set_xlim(0.1, 2.15)\n",
     "ax2.set_ylim(-2.25, 2.25)\n",
-    "ax2.set_title('t='+str(t))\n",
+    "ax2.set_title(\"Nonlinear solve\")\n",
+    "\n",
     "plt.tight_layout()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -870,7 +967,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.17"
+   "version": "3.10.16"
   }
  },
  "nbformat": 4,
diff --git a/freegsnke/jtor_update.py b/freegsnke/jtor_update.py
index b6d9161d..fd4b9e7d 100644
--- a/freegsnke/jtor_update.py
+++ b/freegsnke/jtor_update.py
@@ -853,3 +853,33 @@ def assign_profile_parameter(
             ffp_values_2,
             ffp_sigma,
         ]
+
+
+class GeneralPprimeFFprime(freegs4e.jtor.GeneralPprimeFFprime, Jtor_universal):
+    """FreeGSNKE profile class adapting the FreeGS4E object with the same name,
+    with a few modifications, to:
+    - retain memory of critical point calculation;
+    - deal with limiter plasma configurations
+    """
+
+    def __init__(self, eq, *args, **kwargs):
+        """Instantiates the object.
+
+        Parameters
+        ----------
+        eq : FreeGSNKE Equilibrium object
+            Specifies the domain properties
+        """
+
+        freegs4e.jtor.GeneralPprimeFFprime.__init__(self, *args, **kwargs)
+        Jtor_universal.__init__(self)
+
+        self.profile_parameter = []
+        self.set_masks(eq=eq)
+
+    def assign_profile_parameter(
+        self,
+    ):
+        """Assigns to the profile object new values for the profile parameters"""
+
+        self.profile_parameter = []
diff --git a/freegsnke/linear_solve.py b/freegsnke/linear_solve.py
index 1287306b..8399c68b 100644
--- a/freegsnke/linear_solve.py
+++ b/freegsnke/linear_solve.py
@@ -40,7 +40,6 @@ def __init__(
         Pm1,
         Rm1,
         Mey,
-        # limiter_handler,
         plasma_norm_factor,
         plasma_resistance_1d,
         max_internal_timestep=0.0001,
@@ -137,7 +136,7 @@ def reset_plasma_resistivity(self, plasma_resistance_1d):
             Vector of 2pi resistivity R/dA for all domain grid points in the reduced plasma domain
         """
         self.plasma_resistance_1d = plasma_resistance_1d
-        self.set_linearization_point(None, None, None)
+        self.set_linearization_point(None, None, None, None)
 
     def reset_timesteps(self, max_internal_timestep, full_timestep):
         """Resets the integration timesteps, calling self.solver.set_timesteps
@@ -155,7 +154,7 @@ def reset_timesteps(self, max_internal_timestep, full_timestep):
             full_timestep=full_timestep, max_internal_timestep=max_internal_timestep
         )
 
-    def set_linearization_point(self, dIydI, hatIy0, Myy_hatIy0):
+    def set_linearization_point(self, dIydI, dIydtheta, hatIy0, Myy_hatIy0):
         """Initialises an implicit-Euler solver with the appropriate matrices for the linearised dynamic problem.
 
         Parameters
@@ -165,6 +164,9 @@ def set_linearization_point(self, dIydI, hatIy0, Myy_hatIy0):
             (active coil currents, vessel normal modes, total plasma current divided by plasma_norm_factor).
             These would typically come from having solved the forward Grad-Shafranov problem. Finite difference Jacobian.
             Calculated by the nl_solver object
+        dIydtheta : np.array
+            partial derivatives of plasma-cell currents on the reduced plasma domain with respect to all plasma current density
+            profile parameters
         hatIy0 : np.array
             Plasma current distribution on the reduced plasma domain (1d) of the equilibrium around which the dynamics is linearised.
             This is normalised by the total plasma current, so that this vector sums to 1.
@@ -174,6 +176,8 @@ def set_linearization_point(self, dIydI, hatIy0, Myy_hatIy0):
         """
         if dIydI is not None:
             self.dIydI = dIydI
+        if dIydtheta is not None:
+            self.dIydtheta = dIydtheta
         if hatIy0 is not None:
             self.hatIy0 = hatIy0
         if Myy_hatIy0 is not None:
@@ -207,9 +211,6 @@ def build_Mmatrix(
         where A = dIy/dId
               B = dIy/dIp
               C = dIy/plasmapars
-        Here we take C=0, that is the linearised dynamics does not account for evolving
-        plasma parameters atm.
-
 
         Parameters
         ----------
@@ -253,13 +254,19 @@ def build_Mmatrix(
 
         # build necessary terms to incorporate forcing term from variations of the profile parameters
         # MIdot + RI = V - self.Vm1Rm12Mey_plus@self.dIydpars@d_profiles_pars_dt
-        # if self.dIydpars is not None:
-        #     Pm1Rm1Mey_plus = np.concatenate(
-        #         (self.Pm1Rm1Mey, JMyy[np.newaxis] / nRp), axis=0
-        #     )
-        #     self.forcing_pars_matrix = np.matmul(Pm1Rm1Mey_plus, self.dIydpars)
-
-    def stepper(self, It, active_voltage_vec, d_profiles_pars_dt=None):
+        self.forcing_pars_matrix = None
+        if self.dIydtheta is not None:
+            Pm1Rm1Mey_plus = np.concatenate(
+                (self.Pm1Rm1Mey, self.Myy_hatIy0[np.newaxis] / nRp), axis=0
+            )
+            self.forcing_pars_matrix = np.matmul(Pm1Rm1Mey_plus, self.dIydtheta)
+
+    def stepper(
+        self,
+        It,
+        active_voltage_vec,
+        dtheta_dt,
+    ):
         """Executes the time advancement. Uses the implicit_euler instance.
 
         Parameters
@@ -269,19 +276,21 @@ def stepper(self, It, active_voltage_vec, d_profiles_pars_dt=None):
             (active currents, vessel normal modes, total plasma current divided by normalisation factor)
         active_voltage_vec : np.array
             voltages applied to the active coils
-        d_profiles_pars_dt : np.array
-            time derivative of the profile parameters, not used atm
-        other parameters are passed in as object attributes
+        dtheta_dt : np.array
+            Vector of plasma current density profile parameters derivateives with respect to t.
         """
+
+        # baseline forcing term (from the active coil voltages)
         self.empty_U[: self.n_active_coils] = active_voltage_vec
         self.forcing[:-1] = np.dot(self.Pm1Rm1, self.empty_U)
         self.forcing[-1] = 0.0
 
-        # add forcing term from time derivative of profile parameters
-        if d_profiles_pars_dt is not None:
-            self.forcing -= np.dot(self.forcing_pars_matrix, d_profiles_pars_dt)
+        # additional forcing due to the time derivative of profile parameters
+        if self.forcing_pars_matrix is not None:
+            self.forcing -= np.dot(self.forcing_pars_matrix, dtheta_dt)
 
         Itpdt = self.solver.full_stepper(It, self.forcing)
+
         return Itpdt
 
     def calculate_linear_growth_rate(
diff --git a/freegsnke/nonlinear_solve.py b/freegsnke/nonlinear_solve.py
index 2639ec2d..1894fe6a 100644
--- a/freegsnke/nonlinear_solve.py
+++ b/freegsnke/nonlinear_solve.py
@@ -63,6 +63,7 @@ def __init__(
         mode_removal=True,
         linearize=True,
         dIydI=None,
+        dIydtheta=None,
         target_relative_tolerance_linearization=1e-8,
         target_dIy=1e-3,
         force_core_mask_linearization=False,
@@ -126,6 +127,28 @@ def __init__(
             Threshold value used to include/exclude vessel normal modes.
             Only modes with smaller characteristic frequencies (larger timescales) are retained.
             If None, max_mode_frequency is set based on the input timestep: max_mode_frequency = 1/(5*full_timestep)
+        dIydI : np.array of size (np.sum(plasma_domain_mask), n_metal_modes+1), optional
+            dIydI_(i,j) = d(Iy_i)/d(I_j)
+            This is the jacobian of the plasma current distribution Iy with respect to all
+            independent metal currents (both active and vessel modes) and to the total plasma current
+            This is provided if known, otherwise calculated here at the linearization eq
+        dIydtheta : np.array of size (np.sum(plasma_domain_mask), n_profile_parameters, optional
+            dIydtheta(i,j) = d(Iy_i)/d(theta_j)
+            This is the jacobian of the plasma current distribution Iy with respect to all
+            plasma current desnity profile parameters
+            This is provided if known, otherwise calculated here at the linearization eq
+        automatic_timestep : (float, float) or False, optional, by default False
+            If not False, this overrides inputs full_timestep and max_internal_timestep:
+            the timescales of the linearised problem are used to set the size of the timestep.
+            The input eq and profiles are used to calculate the fastest growthrate, t_growthrate, henceforth,
+            full_timestep = automatic_timestep[0]*t_growthrate
+            max_internal_timestep = automatic_timestep[1]*full_timestep
+        mode_removal : bool, optional, by default True
+            It True, vessel normal modes are dropped after dIydI is calculated
+            Modes that couple with the plasma less than min_dIy_dI than the strongest mode, are dropped.
+            This criterion is applied based on the actual dIydI, calculated on GS solutions.
+        linearize : bool, optional, by default True
+            Whether to set up the linearization of the evolutive problem
         fix_n_vessel_modes : int
             If -1, modes are selected based on max_mode_frequency, threshold_dIy_dI and min_dIy_dI.
             If a non-negative integer, the number of vessel modes is fixed accordingly. max_mode_frequency, threshold_dIy_dI and min_dIy_dI are not used.
@@ -233,6 +256,7 @@ def __init__(
         self.dt_step = full_timestep
         self.max_internal_timestep = max_internal_timestep
         self.set_plasma_resistivity(plasma_resistivity)
+        self.target_dIy = target_dIy
 
         # prepare for mode selection
         if max_mode_frequency is None:
@@ -408,6 +432,48 @@ def __init__(
             (self.approved_target_dIy, [target_dIy])
         )
 
+        # starting dtheta values for Jacobian calculation
+        self.approved_target_dtheta = target_dIy * np.ones(self.n_profiles_parameters)
+        self.starting_dtheta = np.zeros(self.n_profiles_parameters)
+        self.final_dtheta_record = np.zeros(self.n_profiles_parameters)
+
+        # for setting profile parameters starting shifts (for Jacobians)
+        if self.profiles_param is not None:
+            # alpha_m
+            self.starting_dtheta[0] = max(profiles.alpha_m * 1e-2, 1e-2)
+            # alpha_n
+            self.starting_dtheta[1] = max(profiles.alpha_n * 1e-2, 1e-2)
+            # paxis/betap/Beta0
+            if self.profiles_param == "paxis":
+                self.starting_dtheta[2] = max(
+                    getattr(profiles, self.profiles_param) * 1e-2, 1e1
+                )
+            elif self.profiles_param == "betap":
+                self.starting_dtheta[2] = max(
+                    getattr(profiles, self.profiles_param) * 1e-2, 1e-4
+                )
+            elif self.profiles_param == "Beta0":
+                self.starting_dtheta[2] = max(
+                    getattr(profiles, self.profiles_param) * 1e-2, 1e-4
+                )
+
+        else:  # lao
+            # alpha coeffs
+            self.starting_dtheta[0 : self.n_profiles_parameters_alpha] = (
+                self.profiles_parameters_vec[self.profiles_alpha_indices] * 1e-3
+            )
+            self.starting_dtheta[0 : self.n_profiles_parameters_alpha][
+                self.starting_dtheta[0 : self.n_profiles_parameters_alpha] == 0
+            ] = 1e-3
+
+            # beta coeffs
+            self.starting_dtheta[self.n_profiles_parameters_alpha :] = (
+                self.profiles_parameters_vec[self.profiles_beta_indices] * 1e-3
+            )
+            self.starting_dtheta[self.n_profiles_parameters_alpha :][
+                self.starting_dtheta[self.n_profiles_parameters_alpha :] == 0
+            ] = 1e-3
+
         # This solves the system of circuit eqs based on an assumption
         # for the direction of the plasma current distribution at time t+dt
         self.simplified_solver_J1 = simplified_solver_J1(
@@ -479,6 +545,11 @@ def __init__(
         self.dIydI_ICs = dIydI
         self.dIydI = dIydI
 
+        # self.dIydtheta is the Jacobian of the plasma current distribution
+        # with respect to the plasma current density profile parameters
+        self.dIydtheta_ICs = dIydtheta
+        self.dIydtheta = dIydtheta
+
         # initialize and set up the linearization
         # input value for dIydI is used when available
         if automatic_timestep == False:
@@ -497,6 +568,7 @@ def __init__(
                 profiles,
                 target_relative_tolerance_linearization=target_relative_tolerance_linearization,
                 dIydI=dIydI,
+                dIydtheta=dIydtheta,
                 verbose=verbose,
                 force_core_mask_linearization=force_core_mask_linearization,
             )
@@ -672,15 +744,15 @@ def build_dIydI_noGS(
             else:
                 starting_dI[j] = 1.0 * self.final_dI_record[j]
 
-            if verbose:
-                print("mode", j, "; approved starting_dI=", starting_dI[j])
-                print(
-                    "approved relative dIy change=",
-                    rel_ndIy,
-                    "; core_check=",
-                    core_check,
-                )
-                print(" ")
+            # if verbose:
+            #     print("mode", j, "; approved starting_dI=", starting_dI[j])
+            #     print(
+            #         "approved relative dIy change=",
+            #         rel_ndIy,
+            #         "; core_check=",
+            #         core_check,
+            #     )
+            #     print(" ")
             self.dIydI_noGS[:, j] = dIydInoGS
             self.rel_ndIy[j] = rel_ndIy
             # self.final_dI_record[j] = starting_dI[j] * self.accepted_target_dIy[j] / rel_ndIy
@@ -738,12 +810,19 @@ def remove_modes(self, eq, selected_modes_mask):
         )
 
         self.linearised_sol.set_linearization_point(
-            dIydI=self.dIydI, hatIy0=self.blended_hatIy, Myy_hatIy0=self.Myy_hatIy0
+            dIydI=self.dIydI,
+            dIydtheta=self.dIydtheta,
+            hatIy0=self.blended_hatIy,
+            Myy_hatIy0=self.Myy_hatIy0,
         )
 
         self.build_current_vec(self.eq1, self.profiles1)
 
-    def set_linear_solution(self, active_voltage_vec, d_profiles_pars_dt=None):
+    def set_linear_solution(
+        self,
+        active_voltage_vec,
+        dtheta_dt,
+    ):
         """Uses the solver of the linearised problem to set up an initial guess for the nonlinear solver
         for the currents at time t+dt. Uses self.currents_vec as I(t).
         Solves GS at time t+dt using the currents derived from the linearised dynamics.
@@ -752,147 +831,424 @@ def set_linear_solution(self, active_voltage_vec, d_profiles_pars_dt=None):
         ----------
         active_voltage_vec : np.array
             Vector of external voltage applied to the active coils during the timestep.
-        d_profiles_pars_dt : None
-            Does not currently use d_profiles_pars_dt
-            The evolution of the profiles coefficient is not accounted by the linearised dynamical evolution.
+        dtheta_dt : np.array
+            Vector of plasma current density profile parameters derivateives with respect to t.
         """
 
         self.trial_currents = self.linearised_sol.stepper(
             It=self.currents_vec,
             active_voltage_vec=active_voltage_vec,
-            d_profiles_pars_dt=d_profiles_pars_dt,
+            dtheta_dt=dtheta_dt,
         )
         self.assign_currents_solve_GS(self.trial_currents, self.rtol_NK)
         self.trial_plasma_psi = np.copy(self.eq2.plasma_psi)
 
-    # # not used atm, used when building the Jacobian of the plasma current distribution with respect to the coefficients of the profiles object
-    # def prepare_build_dIydpars(self, profiles, rtol_NK, target_dIy, starting_dpars):
-    #     """Prepares to compute the term d(Iy)/d(alpha_m, alpha_n, profifile_par)
-    #     where profiles_par = paxis or betap.
-    #     It infers the value of delta(indep_variable) corresponding to a change delta(I_y)
-    #     with norm(delta(I_y))=target_dIy.
+    def prepare_build_dIydtheta(
+        self,
+        profiles,
+        rtol_NK,
+        target_dIy,
+        starting_dtheta,
+        verbose=False,
+    ):
+        """
 
-    #     Parameters
-    #     ----------
-    #     profiles : FreeGS4E profiles object
-    #         The profiles object of the initial condition equilibrium, i.e. the linearization point.
-    #     rtol_NK : float
-    #         Relative tolerance to be used in the static GS problems.
-    #     target_dIy : float
-    #         Target value for the norm of delta(I_y), on which th finite difference derivative is calculated.
-    #     starting_dpars : tuple (d_alpha_m, d_alpha_n, relative_d_profiles_par)
-    #         Initial value to be used as delta(indep_variable) to infer the slope of norm(delta(I_y))/delta(indep_variable).
-    #         Note that the first two values in the tuple are absolute deltas,
-    #         while the third value is relative, d_profiles_par = relative_d_profiles_par * profiles_par
-    #     """
+        Prepares to compute the term d(Iy)/d(theta) where theta is a vector of the parameters
+        that define the plasma current density profile being used.
+        It infers the value of delta(indep_variable) corresponding to a change delta(I_y)
+        with norm(delta(I_y))=target_dIy.
 
-    #     current_ = np.copy(self.currents_vec)
+        Parameters
+        ----------
+        profiles : FreeGS4E profiles object
+            The profiles object of the initial condition equilibrium, i.e. the linearization point.
+        rtol_NK : float
+            Relative tolerance to be used in the static GS problems.
+        target_dIy : array
+            Target values for the norm of delta(I_y), on which the finite difference derivative is calculated.
+        starting_dtheta : array
+            Initial value to be used as delta(indep_variable) to infer the slope of norm(delta(I_y))/delta(indep_variable).
+        verbose : bool
+            Print some intermediate values, default is Flase.
+        """
 
-    #     # vary alpha_m
-    #     self.check_and_change_profiles(
-    #         profiles_coefficients=(
-    #             profiles.alpha_m + starting_dpars[0],
-    #             profiles.alpha_n,
-    #         )
-    #     )
-    #     self.assign_currents_solve_GS(current_, rtol_NK)
-    #     dIy_0 = self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
-    #     rel_ndIy_0 = np.linalg.norm(dIy_0) / self.nIy
-    #     self.final_dpars_record[0] = starting_dpars[0] * target_dIy / rel_ndIy_0
-
-    #     # vary alpha_n
-    #     self.check_and_change_profiles(
-    #         profiles_coefficients=(
-    #             profiles.alpha_m,
-    #             profiles.alpha_n + starting_dpars[1],
-    #         )
-    #     )
-    #     self.assign_currents_solve_GS(current_, rtol_NK)
-    #     dIy_0 = self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
-    #     rel_ndIy_0 = np.linalg.norm(dIy_0) / self.nIy
-    #     self.final_dpars_record[1] = starting_dpars[1] * target_dIy / rel_ndIy_0
-
-    #     # vary paxis or betap
-    #     self.check_and_change_profiles(
-    #         profiles_coefficients=(profiles.alpha_m, profiles.alpha_n),
-    #         profiles_parameter=(1 + starting_dpars[2]) * profiles.profiles_parameter,
-    #     )
-    #     self.assign_currents_solve_GS(current_, rtol_NK)
-    #     dIy_0 = self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
-    #     rel_ndIy_0 = np.linalg.norm(dIy_0) / self.nIy
-    #     self.final_dpars_record[2] = (
-    #         starting_dpars[2] * profiles.profiles_parameter * target_dIy / rel_ndIy_0
-    #     )
+        current_ = np.copy(self.currents_vec)
 
-    # # not used atm, builds the Jacobian of the plasma current distribution with respect to the coefficients of the profiles object
-    # # can only handle ConstrainPaxisIp and ConstrainBetapIp profiles families.
-    # def build_dIydIpars(self, profiles, rtol_NK, verbose=False):
-    #     """Compute the matrix d(Iy)/d(alpha_m, alpha_n, profifile_par) as a finite difference derivative,
-    #     using the value of delta(indep_viriable) inferred earlier by self.prepare_build_dIypars.
+        # storage
+        dIy_0 = np.zeros((len(self.Iy), self.n_profiles_parameters))
+        rel_ndIy_0 = np.zeros(self.n_profiles_parameters)
 
-    #     Parameters
-    #     ----------
-    #     profiles : FreeGS4E profiles object
-    #         The profiles object of the initial condition equilibrium, i.e. the linearization point.
-    #     rtol_NK : float
-    #         Relative tolerance to be used in the static GS problems.
+        # carry out the initial perturbations
+        if self.profiles_param is not None:
 
-    #     """
+            # vary alpha_m
+            self.check_and_change_profiles(
+                profiles_parameters={
+                    "alpha_m": profiles.alpha_m + starting_dtheta[0],
+                    "alpha_n": profiles.alpha_n,
+                    self.profiles_param: getattr(profiles, self.profiles_param),
+                }
+            )
 
-    #     current_ = np.copy(self.currents_vec)
+            # reset plasma flux map to original
+            self.eq2.plasma_psi = np.copy(self.eq1.plasma_psi)
+            self.assign_currents_solve_GS(current_, rtol_NK)
+            dIy_0[:, 0] = (
+                self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
+            )
+            rel_ndIy_0[0] = np.linalg.norm(dIy_0[:, 0]) / self.nIy
+            self.final_dtheta_record[0] = (
+                starting_dtheta[0] * target_dIy[0] / rel_ndIy_0[0]
+            )
 
-    #     # vary alpha_m
-    #     self.check_and_change_profiles(
-    #         profiles_coefficients=(
-    #             profiles.alpha_m + self.final_dpars_record[0],
-    #             profiles.alpha_n,
-    #         ),
-    #         profiles_parameter=profiles.profiles_parameter,
-    #     )
-    #     self.assign_currents_solve_GS(current_, rtol_NK)
-    #     dIy_1 = self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
-    #     self.dIydpars[:, 0] = dIy_1 / self.final_dpars_record[0]
-    #     if verbose:
-    #         print(
-    #             "alpha_m gradient calculated on the finite difference: delta_alpha_m =",
-    #             self.final_dpars_record[0],
-    #             ", norm(deltaIy) =",
-    #             np.linalg.norm(dIy_1),
-    #         )
-
-    #     # vary alpha_n
-    #     self.check_and_change_profiles(
-    #         profiles_coefficients=(
-    #             profiles.alpha_m,
-    #             profiles.alpha_n + self.final_dpars_record[1],
-    #         )
-    #     )
-    #     self.assign_currents_solve_GS(current_, rtol_NK)
-    #     dIy_1 = self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
-    #     self.dIydpars[:, 1] = dIy_1 / self.final_dpars_record[1]
-    #     if verbose:
-    #         print(
-    #             "alpha_n gradient calculated on the finite difference: delta_alpha_n =",
-    #             self.final_dpars_record[1],
-    #             ", norm(deltaIy) =",
-    #             np.linalg.norm(dIy_1),
-    #         )
-
-    #     # vary paxis or betap
-    #     self.check_and_change_profiles(
-    #         profiles_coefficients=(profiles.alpha_m, profiles.alpha_n),
-    #         profiles_parameter=profiles.profiles_parameter + self.final_dpars_record[2],
-    #     )
-    #     self.assign_currents_solve_GS(current_, rtol_NK)
-    #     dIy_1 = self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
-    #     self.dIydpars[:, 2] = dIy_1 / self.final_dpars_record[2]
-    #     if verbose:
-    #         print(
-    #             "profiles_par gradient calculated on the finite difference: delta_profiles_par =",
-    #             self.final_dpars_record[2],
-    #             ", norm(deltaIy) =",
-    #             np.linalg.norm(dIy_1),
-    #         )
+            if verbose:
+                print("")
+                print("Profile parameter: alpha_m:")
+                print(f"  Initial delta parameter = {starting_dtheta[0]}")
+                print(f"  Initial relative Iy change = {rel_ndIy_0[0]}")
+                print(f"  Final delta parameter = {self.final_dtheta_record[0]}")
+
+            # vary alpha_n
+            self.check_and_change_profiles(
+                profiles_parameters={
+                    "alpha_m": profiles.alpha_m,
+                    "alpha_n": profiles.alpha_n + starting_dtheta[1],
+                    self.profiles_param: getattr(profiles, self.profiles_param),
+                }
+            )
+
+            # reset plasma flux map to original
+            self.eq2.plasma_psi = np.copy(self.eq1.plasma_psi)
+            self.assign_currents_solve_GS(current_, rtol_NK)
+            dIy_0[:, 1] = (
+                self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
+            )
+            rel_ndIy_0[1] = np.linalg.norm(dIy_0[:, 1]) / self.nIy
+            self.final_dtheta_record[1] = (
+                starting_dtheta[1] * target_dIy[1] / rel_ndIy_0[1]
+            )
+
+            if verbose:
+                print("")
+                print("Profile parameter: alpha_n:")
+                print(f"  Initial delta parameter = {starting_dtheta[1]}")
+                print(f"  Initial relative Iy change = {rel_ndIy_0[1]}")
+                print(f"  Final delta parameter = {self.final_dtheta_record[1]}")
+
+            # vary paxis, betap or Beta0
+            self.check_and_change_profiles(
+                profiles_parameters={
+                    "alpha_m": profiles.alpha_m,
+                    "alpha_n": profiles.alpha_n,
+                    self.profiles_param: getattr(profiles, self.profiles_param)
+                    + starting_dtheta[2],
+                }
+            )
+
+            # reset plasma flux map to original
+            self.eq2.plasma_psi = np.copy(self.eq1.plasma_psi)
+            self.assign_currents_solve_GS(current_, rtol_NK)
+            dIy_0[:, 2] = (
+                self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
+            )
+            rel_ndIy_0[2] = np.linalg.norm(dIy_0[:, 2]) / self.nIy
+            self.final_dtheta_record[2] = (
+                starting_dtheta[2] * target_dIy[2] / rel_ndIy_0[2]
+            )
+
+            if verbose:
+                print("")
+                print(f"Profile parameter: {self.profiles_param}:")
+                print(f"  Initial delta parameter = {starting_dtheta[2]}")
+                print(f"  Initial relative Iy change = {rel_ndIy_0[2]}")
+                print(f"  Final delta parameter = {self.final_dtheta_record[2]}")
+
+            # reset profiles in profiles1 and profiles2 objects
+            self.check_and_change_profiles(
+                profiles_parameters={
+                    "alpha_m": profiles.alpha_m,
+                    "alpha_n": profiles.alpha_n,
+                    self.profiles_param: getattr(profiles, self.profiles_param),
+                }
+            )
+
+        else:  # this is particular to the Lao profile coefficients (which there may be few or many of)
+
+            # for each alpha coefficient
+            alpha_base = profiles.alpha.copy()
+            for i in range(0, self.n_profiles_parameters_alpha):
+                alpha_shift = alpha_base.copy()
+                alpha_shift[i] += starting_dtheta[i]  # perturb the term
+                if profiles.alpha_logic:
+                    alpha_shift[-1] -= starting_dtheta[
+                        i
+                    ]  # final alpha term needs perturbing too
+
+                self.check_and_change_profiles(
+                    profiles_parameters={
+                        "alpha": alpha_shift,
+                        "beta": profiles.beta,
+                    }
+                )
+
+                # reset plasma flux map to original
+                self.eq2.plasma_psi = np.copy(self.eq1.plasma_psi)
+                self.assign_currents_solve_GS(current_, rtol_NK)
+                dIy_0[:, i] = (
+                    self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
+                )
+                rel_ndIy_0[i] = np.linalg.norm(dIy_0[:, i]) / self.nIy
+                self.final_dtheta_record[i] = (
+                    starting_dtheta[i] * target_dIy[i] / rel_ndIy_0[i]
+                )
+
+                if verbose:
+                    print("")
+                    print(f"Profile parameter: alpha_{i}:")
+                    print(f"  Initial delta parameter = {starting_dtheta[i]}")
+                    print(f"  Initial relative Iy change = {rel_ndIy_0[i]}")
+                    print(f"  Final delta parameter = {self.final_dtheta_record[i]}")
+
+            # for each beta coefficient
+            beta_base = profiles.beta.copy()
+            for i in range(0, self.n_profiles_parameters_beta):
+                beta_shift = beta_base.copy()
+                beta_shift[i] += starting_dtheta[
+                    i + self.n_profiles_parameters_alpha
+                ]  # perturb the term required
+                if profiles.beta_logic:
+                    beta_shift[-1] -= starting_dtheta[
+                        i + self.n_profiles_parameters_alpha
+                    ]  # final beta term needs perturbing too
+                self.check_and_change_profiles(
+                    profiles_parameters={
+                        "alpha": profiles.alpha,
+                        "beta": beta_shift,
+                    }
+                )
+
+                # reset plasma flux map to original
+                self.eq2.plasma_psi = np.copy(self.eq1.plasma_psi)
+                self.assign_currents_solve_GS(current_, rtol_NK)
+                dIy_0[:, i + self.n_profiles_parameters_alpha] = (
+                    self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
+                )
+                rel_ndIy_0[i + self.n_profiles_parameters_alpha] = (
+                    np.linalg.norm(dIy_0[:, i + self.n_profiles_parameters_alpha])
+                    / self.nIy
+                )
+                self.final_dtheta_record[i + self.n_profiles_parameters_alpha] = (
+                    starting_dtheta[i + self.n_profiles_parameters_alpha]
+                    * target_dIy[i]
+                    / rel_ndIy_0[i + self.n_profiles_parameters_alpha]
+                )
+
+                if verbose:
+                    print("")
+                    print(f"Profile parameter: beta_{i}:")
+                    print(
+                        f"  Initial delta parameter = {starting_dtheta[i+self.n_profiles_parameters_alpha]}"
+                    )
+                    print(
+                        f"  Initial relative Iy change = {rel_ndIy_0[i+self.n_profiles_parameters_alpha]}"
+                    )
+                    print(
+                        f"  Final delta parameter = {self.final_dtheta_record[i+self.n_profiles_parameters_alpha]}"
+                    )
+
+                # reset profiles in profiles1 and profiles2 objects
+                self.check_and_change_profiles(
+                    profiles_parameters={
+                        "alpha": profiles.alpha,
+                        "beta": profiles.beta,
+                    }
+                )
+
+        return dIy_0 / starting_dtheta, rel_ndIy_0
+
+    def build_dIydtheta(self, profiles, rtol_NK, verbose=False):
+        """
+
+        Compute the matrix d(Iy)/d(theta) as a finite difference derivative,
+        using the value of delta(indep_viriable) inferred earlier by self.prepare_build_dIypars.
+
+        Parameters
+        ----------
+        profiles : FreeGS4E profiles object
+            The profiles object of the initial condition equilibrium, i.e. the linearization point.
+        rtol_NK : float
+            Relative tolerance to be used in the static GS problems.
+
+        """
+
+        # the final perturbations to calc the Jacobian with
+        final_theta = 1.0 * self.final_dtheta_record
+
+        current_ = np.copy(self.currents_vec)
+
+        # storage
+        dIydtheta = np.zeros((len(self.Iy), self.n_profiles_parameters))
+        rel_ndIy = np.zeros(self.n_profiles_parameters)
+
+        # carry out the initial perturbations
+        if self.profiles_param is not None:
+
+            # vary alpha_m
+            self.check_and_change_profiles(
+                profiles_parameters={
+                    "alpha_m": profiles.alpha_m + final_theta[0],
+                    "alpha_n": profiles.alpha_n,
+                    self.profiles_param: getattr(profiles, self.profiles_param),
+                }
+            )
+            # reset plasma flux map to original
+            self.eq2.plasma_psi = np.copy(self.eq1.plasma_psi)
+            self.assign_currents_solve_GS(current_, rtol_NK)
+            dIy_1 = self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
+            rel_ndIy[0] = np.linalg.norm(dIy_1) / self.nIy
+            dIydtheta[:, 0] = dIy_1 / final_theta[0]
+            if verbose:
+                print("")
+                print(f"Profile parameter: alpha_m:")
+                print(f"  Final relative Iy change = {rel_ndIy[0]}")
+                print(
+                    f"  Initial vs. Final GS residual: {self.NK.initial_rel_residual} vs. {self.NK.relative_change}"
+                )
+
+            # vary alpha_n
+            self.check_and_change_profiles(
+                profiles_parameters={
+                    "alpha_m": profiles.alpha_m,
+                    "alpha_n": profiles.alpha_n + final_theta[1],
+                    self.profiles_param: getattr(profiles, self.profiles_param),
+                }
+            )
+            # reset plasma flux map to original
+            self.eq2.plasma_psi = np.copy(self.eq1.plasma_psi)
+            self.assign_currents_solve_GS(current_, rtol_NK)
+            dIy_1 = self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
+            rel_ndIy[1] = np.linalg.norm(dIy_1) / self.nIy
+            dIydtheta[:, 1] = dIy_1 / final_theta[1]
+            if verbose:
+                print("")
+                print(f"Profile parameter: alpha_n:")
+                print(f"  Final relative Iy change = {rel_ndIy[1]}")
+                print(
+                    f"  Initial vs. Final GS residual: {self.NK.initial_rel_residual} vs. {self.NK.relative_change}"
+                )
+
+            # vary paxis, betap or Beta0
+            self.check_and_change_profiles(
+                profiles_parameters={
+                    "alpha_m": profiles.alpha_m,
+                    "alpha_n": profiles.alpha_n,
+                    self.profiles_param: getattr(profiles, self.profiles_param)
+                    + final_theta[2],
+                }
+            )
+            # reset plasma flux map to original
+            self.eq2.plasma_psi = np.copy(self.eq1.plasma_psi)
+            self.assign_currents_solve_GS(current_, rtol_NK)
+            dIy_1 = self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
+            rel_ndIy[2] = np.linalg.norm(dIy_1) / self.nIy
+            dIydtheta[:, 2] = dIy_1 / final_theta[2]
+            if verbose:
+                print("")
+                print(f"Profile parameter: {self.profiles_param}:")
+                print(f"  Final relative Iy change = {rel_ndIy[2]}")
+                print(
+                    f"  Initial vs. Final GS residual: {self.NK.initial_rel_residual} vs. {self.NK.relative_change}"
+                )
+
+            # reset profiles in profiles1 and profiles2 objects
+            self.check_and_change_profiles(
+                profiles_parameters={
+                    "alpha_m": profiles.alpha_m,
+                    "alpha_n": profiles.alpha_n,
+                    self.profiles_param: getattr(profiles, self.profiles_param),
+                }
+            )
+
+        else:  # this is particular to the Lao profile coefficients (which there may be few or many of)
+
+            # for each alpha coefficient
+
+            alpha_base = profiles.alpha.copy()
+            for i in range(0, self.n_profiles_parameters_alpha):
+                alpha_shift = alpha_base.copy()
+                alpha_shift[i] += final_theta[i]  # perturb the term
+                if profiles.alpha_logic:
+                    alpha_shift[-1] -= final_theta[
+                        i
+                    ]  # final alpha term needs perturbing too
+
+                self.check_and_change_profiles(
+                    profiles_parameters={
+                        "alpha": alpha_shift,
+                        "beta": profiles.beta,
+                    }
+                )
+                # reset plasma flux map to original
+                self.eq2.plasma_psi = np.copy(self.eq1.plasma_psi)
+                self.assign_currents_solve_GS(current_, rtol_NK)
+                dIy_1 = self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
+                rel_ndIy[i] = np.linalg.norm(dIy_1) / self.nIy
+                dIydtheta[:, i] = dIy_1 / final_theta[i]
+                if verbose:
+                    print("")
+                    print(f"Profile parameter: alpha_{i}:")
+                    print(f"  Final relative Iy change = {rel_ndIy[i]}")
+                    print(
+                        f"  Initial vs. Final GS residual: {self.NK.initial_rel_residual} vs. {self.NK.relative_change}"
+                    )
+
+            # for each beta coefficient
+            beta_base = profiles.beta.copy()
+            for i in range(0, self.n_profiles_parameters_beta):
+                beta_shift = beta_base.copy()
+                beta_shift[i] += final_theta[
+                    i + self.n_profiles_parameters_alpha
+                ]  # perturb the term required
+                if profiles.beta_logic:
+                    beta_shift[-1] -= final_theta[
+                        i + self.n_profiles_parameters_alpha
+                    ]  # final beta term needs perturbing too
+                self.check_and_change_profiles(
+                    profiles_parameters={
+                        "alpha": profiles.alpha,
+                        "beta": beta_shift,
+                    }
+                )
+                # reset plasma flux map to original
+                self.eq2.plasma_psi = np.copy(self.eq1.plasma_psi)
+                self.assign_currents_solve_GS(current_, rtol_NK)
+                dIy_1 = self.limiter_handler.Iy_from_jtor(self.profiles2.jtor) - self.Iy
+                rel_ndIy[i + self.n_profiles_parameters_alpha] = (
+                    np.linalg.norm(dIy_1) / self.nIy
+                )
+                dIydtheta[:, i + self.n_profiles_parameters_alpha] = (
+                    dIy_1 / final_theta[i + self.n_profiles_parameters_alpha]
+                )
+
+                if verbose:
+                    print("")
+                    print(f"Profile parameter: beta_{i}:")
+                    print(
+                        f"  Final relative Iy change = {rel_ndIy[i+self.n_profiles_parameters_alpha]}"
+                    )
+                    print(
+                        f"  Initial vs. Final GS residual: {self.NK.initial_rel_residual} vs. {self.NK.relative_change}"
+                    )
+
+            # reset profiles in profiles1 and profiles2 objects
+            self.check_and_change_profiles(
+                profiles_parameters={
+                    "alpha": profiles.alpha,
+                    "beta": profiles.beta,
+                }
+            )
+
+        return dIydtheta, rel_ndIy
 
     def prepare_build_dIydI_j(
         self,
@@ -986,11 +1342,14 @@ def build_linearization(
         eq,
         profiles,
         dIydI,
+        dIydtheta,
         target_relative_tolerance_linearization,
         force_core_mask_linearization,
         verbose,
     ):
-        """Builds the Jacobian d(Iy)/dI to set up the solver of the linearised problem.
+        """
+
+        Builds the Jacobians d(Iy)/dI and d(Iy)/dtheta to set up the solver of the linearised problem.
 
         Parameters
         ----------
@@ -999,17 +1358,19 @@ def build_linearization(
         profiles : FreeGS4E profiles Object
             profiles properties of the equilibrium around which to linearise.
         dIydI : np.array
-            input Jacobian, enter where available, otherwise this will be calculated here
-        rtol_NK : float
+            input Jacobian of plasma current density to metal currents, enter where available, otherwise this will be calculated here
+        dIydtheta : np.array
+            input Jacobian of plasma current density to pasma current density profile parameters, enter where available, otherwise this will be calculated here
+        target_relative_tolerance_linearization : float
             Relative tolerance to be used in the static GS problems.
-        target_dIy : float, by default 0.001.
-            Target relative value for the norm of delta(I_y), on which the finite difference derivative is calculated.
+        force_core_mask_linearization : bool
+            Forces the perturbed plasma core mask to remain the same as the original when calculating the dIydI Jacobian.
         """
 
-        if (dIydI is None) and (self.dIydI is None):
-            self.build_current_vec(eq, profiles)
-            self.Iy = self.limiter_handler.Iy_from_jtor(profiles.jtor).copy()
-            self.nIy = np.linalg.norm(self.Iy)
+        # if (dIydI is None) and (self.dIydI is None):
+        self.build_current_vec(eq, profiles)
+        self.Iy = self.limiter_handler.Iy_from_jtor(profiles.jtor).copy()
+        self.nIy = np.linalg.norm(self.Iy)
 
         self.R0 = eq.Rcurrent()
         self.Z0 = eq.Zcurrent()
@@ -1019,6 +1380,7 @@ def build_linearization(
         # this uses the variation to Jtor caused by the coil's contribution to the flux, ignoring the response of the plasma
 
         # build/update dIydI
+        # dIydI = 1
         if dIydI is None:
             if self.dIydI_ICs is None:
                 print(
@@ -1073,19 +1435,6 @@ def build_linearization(
                                 == 0
                             )
 
-                    if verbose:
-                        print("mode", j)
-                        print(
-                            "Initial delta_current=",
-                            self.starting_dI[j],
-                        )
-                        print(
-                            "Initial relative Iy change=",
-                            ndIy,
-                            "; core_check=",
-                            core_check,
-                        )
-
                     if (
                         np.abs(np.log10(self.final_dI_record[j] / self.starting_dI[j]))
                         > 0.5
@@ -1127,27 +1476,16 @@ def build_linearization(
                         self.final_dI_record[j] = 1.0 * self.starting_dI[j]
 
                     if verbose:
+                        print("")
+                        print(f"Mode: {j}")
+                        print(f"  Initial delta_current = {self.starting_dI[j]}")
+                        print(f"  Initial relative Iy change = {ndIy}")
+                        print(f"  Final delta_current = {self.final_dI_record[j]}")
+                        print("")
+                        print(f"  Final relative Iy change = {rel_ndIy}")
                         print(
-                            "Final current shift=",
-                            self.final_dI_record[j],
+                            f"  Initial vs. Final GS residual: {self.NK.initial_rel_residual} vs. {self.NK.relative_change}"
                         )
-                        # if force_core_mask_linearization:
-                        try:
-                            print(
-                                "Final relative Iy change=",
-                                rel_ndIy,
-                                "; core_check=",
-                                core_check,
-                            )
-                        except:
-                            print("failed printout!")
-                        print(
-                            "Initial residual=",
-                            self.NK.initial_rel_residual,
-                            ". Final residual=",
-                            self.NK.relative_change,
-                        )
-                        print(" ")
 
                     self.dIydI[:, j] = np.copy(dIydIj)
                     self.psideltaI[j] = np.copy(self.eq2.psi())
@@ -1163,6 +1501,55 @@ def build_linearization(
             self.dIydI = dIydI
             self.dIydI_ICs = np.copy(self.dIydI)
 
+        # compose the vector of initial delta_theta (profile parameters) to be used for the finite difference calculation
+        # - this uses the variation to Jtor caused by the parameter's contribution to the change in poloidal flux, ignoring the response of the plasma
+
+        # build/update dIydtheta
+        if dIydtheta is None and force_core_mask_linearization is False:
+            if self.dIydtheta_ICs is None:
+                print(
+                    f"Building the {self.plasma_domain_size} x {self.n_profiles_parameters} Jacobian (dIy/dtheta)",
+                    "of plasma current density (inside the LCFS)",
+                    "with respect to all plasma current density profile parameters within Jtor.",
+                )
+
+                self.dIydtheta = np.zeros(
+                    (self.plasma_domain_size, self.n_profiles_parameters)
+                )
+
+                profiles_copy = deepcopy(profiles)
+
+                # prepare to build the Jacobian by finding appropriate step size
+                dIydtheta, ndIy = self.prepare_build_dIydtheta(
+                    profiles=profiles_copy,
+                    rtol_NK=target_relative_tolerance_linearization,
+                    target_dIy=self.approved_target_dtheta,
+                    starting_dtheta=self.starting_dtheta,
+                    verbose=verbose,
+                )
+
+                if (
+                    np.abs(np.log10(self.final_dtheta_record / self.starting_dtheta))
+                    > 0.5
+                ).any():
+
+                    dIydtheta, rel_ndIy = self.build_dIydtheta(
+                        profiles=profiles_copy,
+                        rtol_NK=target_relative_tolerance_linearization,
+                        verbose=verbose,
+                    )
+
+                else:
+                    self.final_dtheta_record = 1.0 * self.starting_dtheta
+
+                self.dIydtheta = np.copy(dIydtheta)
+                self.dIydtheta_ICs = np.copy(self.dIydtheta)
+            else:
+                self.dIydtheta = np.copy(self.dIydtheta_ICs)
+        else:
+            self.dIydtheta = dIydtheta
+            self.dIydtheta_ICs = np.copy(self.dIydtheta)
+
     def set_plasma_resistivity(self, plasma_resistivity):
         """Function to set the resistivity of the plasma.
         self.plasma_resistance_1d is the diagonal of the matrix R_yy, the plasma resistance matrix.
@@ -1294,25 +1681,66 @@ def get_profiles_values(self, profiles):
         # note the parameters here are the same that should be provided
         # to the stepper if these are time evolving
         if self.profiles_type == "ConstrainPaxisIp":
+            self.n_profiles_parameters = 3
             self.profiles_parameters = {
-                "paxis": profiles.paxis,
                 "alpha_m": profiles.alpha_m,
                 "alpha_n": profiles.alpha_n,
+                "paxis": profiles.paxis,
             }
+            self.profiles_param = "paxis"
+            self.profiles_parameters_vec = np.array(
+                [profiles.alpha_m, profiles.alpha_n, profiles.paxis]
+            )
         elif self.profiles_type == "ConstrainBetapIp":
+            self.n_profiles_parameters = 3
             self.profiles_parameters = {
-                "betap": profiles.betap,
                 "alpha_m": profiles.alpha_m,
                 "alpha_n": profiles.alpha_n,
+                "betap": profiles.betap,
             }
+            self.profiles_param = "betap"
+            self.profiles_parameters_vec = np.array(
+                [profiles.alpha_m, profiles.alpha_n, profiles.betap]
+            )
         elif self.profiles_type == "Fiesta_Topeol":
+            self.n_profiles_parameters = 3
             self.profiles_parameters = {
-                "beta0": profiles.beta0,
                 "alpha_m": profiles.alpha_m,
                 "alpha_n": profiles.alpha_n,
+                "Beta0": profiles.Beta0,
             }
+            self.profiles_param = "Beta0"
+            self.profiles_parameters_vec = np.array(
+                [profiles.alpha_m, profiles.alpha_n, profiles.Beta0]
+            )
         elif self.profiles_type == "Lao85":
+            self.n_profiles_parameters_alpha = len(profiles.alpha)
+            self.n_profiles_parameters_beta = len(profiles.beta)
+            if profiles.alpha_logic:
+                self.n_profiles_parameters_alpha -= 1
+            if profiles.beta_logic:
+                self.n_profiles_parameters_beta -= 1
+            self.n_profiles_parameters = (
+                self.n_profiles_parameters_alpha + self.n_profiles_parameters_beta
+            )
+
+            self.profiles_alpha_indices = slice(0, self.n_profiles_parameters_alpha)
+            alpha_shift = 0
+            if profiles.alpha_logic:
+                alpha_shift += 1
+
+            self.profiles_beta_indices = slice(
+                self.n_profiles_parameters_alpha + alpha_shift,
+                self.n_profiles_parameters_alpha
+                + alpha_shift
+                + self.n_profiles_parameters_beta,
+            )
+
             self.profiles_parameters = {"alpha": profiles.alpha, "beta": profiles.beta}
+            self.profiles_param = None
+            self.profiles_parameters_vec = np.concatenate(
+                (profiles.alpha, profiles.beta)
+            )
 
     def get_vessel_currents(self, eq):
         """Uses the input equilibrium to extract values for all metal currents,
@@ -1325,9 +1753,6 @@ def get_vessel_currents(self, eq):
             Initial equilibrium. eq.tokamak is used to extract current values.
         """
         self.vessel_currents_vec = eq.tokamak.getCurrentsVec()
-        # eq_currents = eq.tokamak.getCurrents()
-        # for i, labeli in enumerate(self.coils_order):
-        #     self.vessel_currents_vec[i] = eq_currents[labeli]
 
     def build_current_vec(self, eq, profiles):
         """Builds the vector of currents in which the dynamics is actually solved, self.currents_vec
@@ -1361,6 +1786,7 @@ def initialize_from_ICs(
         profiles,
         target_relative_tolerance_linearization=1e-7,
         dIydI=None,
+        dIydtheta=None,
         force_core_mask_linearization=False,
         verbose=False,
     ):
@@ -1446,6 +1872,7 @@ def initialize_from_ICs(
             self.eq1,
             self.profiles1,
             dIydI=dIydI,
+            dIydtheta=dIydtheta,
             target_relative_tolerance_linearization=target_relative_tolerance_linearization,
             force_core_mask_linearization=force_core_mask_linearization,
             verbose=verbose,
@@ -1457,7 +1884,10 @@ def initialize_from_ICs(
         # transfer linearization to linear solver:
         self.Myy_hatIy0 = self.handleMyy.dot(self.hatIy)
         self.linearised_sol.set_linearization_point(
-            dIydI=self.dIydI_ICs, hatIy0=self.blended_hatIy, Myy_hatIy0=self.Myy_hatIy0
+            dIydI=self.dIydI_ICs,
+            dIydtheta=self.dIydtheta_ICs,
+            hatIy0=self.blended_hatIy,
+            Myy_hatIy0=self.Myy_hatIy0,
         )
 
     def step_complete_assign(self, working_relative_tol_GS, from_linear=False):
@@ -1861,9 +2291,9 @@ def check_and_change_profiles(self, profiles_parameters=None):
             for par in profiles_parameters:
                 setattr(self.profiles1, par, profiles_parameters[par])
                 setattr(self.profiles2, par, profiles_parameters[par])
-                if self.profiles_type == "Lao85":
-                    self.profiles1.initialize_profile()
-                    self.profiles2.initialize_profile()
+            if self.profiles_type == "Lao85":
+                self.profiles1.initialize_profile()
+                self.profiles2.initialize_profile()
             self.profiles_change_flag = 1
 
     def nlstepper(
@@ -1974,12 +2404,34 @@ def nlstepper(
             NK iterations are interrupted when this limit is surpassed.
         """
 
+        # retrieve the old profile parameter values
+        self.get_profiles_values(self.profiles1)
+        old_params = self.profiles_parameters_vec
+
         # check if profiles parameters are being evolved
         # and action the change where necessary
         self.check_and_change_profiles(
             profiles_parameters=profiles_parameters,
         )
 
+        # retrieve the new profile parameter values (if present)
+        self.get_profiles_values(self.profiles1)
+        new_params = self.profiles_parameters_vec
+
+        # calculate change in profiles across timestep: (profiles(t+dt)-profiles(t))/dt
+        # should be zero if no change
+        if self.profiles_type == "Lao85":
+            old_alphas = old_params[self.profiles_alpha_indices]
+            old_betas = old_params[self.profiles_beta_indices]
+            new_alphas = new_params[self.profiles_alpha_indices]
+            new_betas = new_params[self.profiles_beta_indices]
+            dtheta_dt = (
+                np.concatenate((new_alphas, new_betas))
+                - np.concatenate((old_alphas, old_betas))
+            ) / self.dt_step
+        else:
+            dtheta_dt = (new_params - old_params) / self.dt_step
+
         # check if plasma resistivity is being evolved
         # and action the change where necessary
         self.check_and_change_plasma_resistivity(
@@ -1989,7 +2441,10 @@ def nlstepper(
         # solves the linearised problem for the currents and assigns
         # results in preparation for the nonlinear calculations
         # Solution and GS equilibrium are assigned to self.trial_currents and self.trial_plasma_psi
-        self.set_linear_solution(active_voltage_vec)
+        self.set_linear_solution(
+            active_voltage_vec=active_voltage_vec,
+            dtheta_dt=dtheta_dt,
+        )
 
         # check Matrix is still applicable
         myy_flag = self.handleMyy.check_Myy(self.hatIy)