Is there any way to solve a ODEs system in Unity the way I need?

Question

I want to "grow" a PreFab in unity using some parameters that are the numerical solution of a ODEs system. I have two scripts: one that literally grows the prefab and another one that solves an ODE using 4th order Runge-Kutta numerical method.

The code for ODE solution is this one:

    using System;
    namespace RungeKutta4
    {
        class Program
        {
            static void Main(string[] args)
            {
                //Incrementers to pass into the known solution
                double t = 0.0;
                double T = 1.0;
                double dt = 0.1;

                // Assign the number of elements needed for the arrays
                int n = (int)(((T - t) / dt)) + 1;

                // Initialize the arrays for the time index 's' and estimates 'annualGrowth' at each index 'i'
                double[] annualGrowth = new double[n]; //y
                double[] s = new double[n]; //t

                // RK4 Variables
                double dy1;
                double dy2;
                double dy3;
                double dy4;

                // RK4 Initializations condizioni iniziali
                int i = 0;
                s[i] = 0.0; //t0 = 0
                annualGrowth[i] = 1.0; //y(t0) = y(0) = 1

                // Iterate and implement the Rk4 Algorithm
                while (i < annualGrowth.Length - 1)
                {

                    dy1 = dt * equation(s[i], annualGrowth[i]);
                    dy2 = dt * equation(s[i] + dt / 2, annualGrowth[i] + dy1 / 2);
                    dy3 = dt * equation(s[i] + dt / 2, annualGrowth[i] + dy2 / 2);
                    dy4 = dt * equation(s[i] + dt, annualGrowth[i] + dy3);

                    s[i + 1] = s[i] + dt;
                    annualGrowth[i + 1] = annualGrowth[i] + (dy1 + 2 * dy2 + 2 * dy3 + dy4) / 6;

                    double t_rounded = Math.Round(t + dt, 2);

                    if (t_rounded % 1 == 0)
                    {
                       Console.WriteLine(" y(" + t_rounded + ")" + " " .PadRight(10) + annualGrowth[i + 1]);
                    }

                    i++;
                    t += dt; //t = t + dt


                };//End Rk4                
                Console.ReadLine();
            }

            // Differential Equation
            public static double equation(double t, double annualGrowth)
            {
                double y_prime;
                double k = 0.1;
                double maxHeight = 5.0;
                return y_prime = k * annualGrowth * (1-(annualGrowth/maxHeight));
            }
        }
    }

Then I have the "growing" script in unity:

    public class Growth : MonoBehaviour
{
    SetupScene setup;
    TimeManager timeManager;
    float currentLenght;
    float currentWidth;
    Genetic geneticInfo;
    bool isMaxLenght; 
    bool isMaxWeigth;
    GameObject newInternode;
    float annualLengthGrowth = 0.05f; 
    float annualWidthGrowth = 0.01f;
    Vector3 annualGrowth;
    bool hasChild;


    // Start is called before the first frame update
    void Awake()
    {
        annualGrowth = new Vector3(annualWidthGrowth, annualLengthGrowth, annualWidthGrowth);
        setup = GameObject.Find("Setup").GetComponent<SetupScene>();
        timeManager = GameObject.Find("Time Manager").GetComponent<TimeManager>();
        currentLenght = 0f;
        currentWidth = 0f;
        geneticInfo = new Genetic(1);
    }

    // Update is called once per frame
    void Update()
    {
        Debug.Log("ciao, sono nato " + gameObject.name);
        Debug.Log(timeManager.IsYearElapsed + " " + gameObject.name);
        if (timeManager.IsYearElapsed)
        {

            float scaleFactorXZ = 1f;
            float scaleFactorY = 1f;
            if (gameObject.name != "Seed")
            {
                scaleFactorXZ = geneticInfo.InternodeMaxWidtht;
                scaleFactorY = geneticInfo.InternodeMaxLenght;
            }
            currentLenght = (gameObject.transform.lossyScale.y + annualGrowth.y);
            currentWidth = (gameObject.transform.lossyScale.x + annualGrowth.x);



            if(currentWidth > geneticInfo.InternodeMaxWidtht)
            {
                currentWidth = geneticInfo.InternodeMaxWidtht;
                isMaxWeigth = true;
            }

            if(currentLenght > geneticInfo.InternodeMaxLenght)
            {
                currentLenght = geneticInfo.InternodeMaxLenght;
                isMaxLenght = true;
            }



            gameObject.transform.localScale = new Vector3(currentWidth/scaleFactorXZ, currentLenght/scaleFactorY, currentWidth/scaleFactorXZ);



            if (isMaxLenght)
            {

                if (!hasChild && setup.Manager.InternodeNumber < geneticInfo.MaxInternodeNumber)
                {
                    setup.Manager.InternodeNumber++;
                    Debug.Log("nuovo internodo creato");
                    setup.LastInternodePosition = gameObject.transform.position + new Vector3(0f, geneticInfo.InternodeMaxLenght, 0f);
                    newInternode = Instantiate(setup.Prefab, new Vector3(setup.LastInternodePosition.x, setup.LastInternodePosition.y, setup.LastInternodePosition.z), Quaternion.identity);
                    newInternode.transform.localScale = Vector3.zero;
                    newInternode.transform.parent = setup.LastInternode.transform;
                    setup.LastInternode = newInternode;

                    newInternode.AddComponent<Growth>();
                    hasChild = true;
                }

                if (isMaxWeigth)
                {
                    gameObject.GetComponent<Growth>().enabled = false;
                }

            }
        }
    }
}

I want to link the numerical solution of ODE script to the value float annualLengthGrowth = 0.05f;in the unity script that is, for now, hardcoded. The problem is that at each update (or timestep) the initial value of ODE system must be the final numerical calculated value of the previous step.

For example, given that at first timestep, the initial value of annualGrowth in ODE script is 1.0, the final calculated value is (for example) annualGrowthT1 = 2.35f, I have to use this value in Unity script at first time step as annualLengthGrowth. At next time step, the initial value for ODE script should be annualGrowthT1 = 2.35f, then the calculated final value should be (for example) annualGrowthT2 = 3.56f and this value must be used in the second time step in unity script as new annualLengthGrowth. And so on.

It should be clear thank to this picture:

Answer 1

I tried to put all together but the unity console gives me the error:

error CS0161: 'Growth.Runge.runge(float, float, float, float, Growth.Runge.Function)': not all code paths return a value.

I thin I missed a return somewhere in the code but I don't know where. Here it is my new code:

using System;
using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class Growth : MonoBehaviour
{
    SetupScene setup; //mi serve per prendere il prefab
    TimeManager timeManager;
    float currentLenght;
    float currentWidth;
    Genetic geneticInfo;
    bool isMaxLenght;  //internodo ha raggiunto la lunghezza massima
    bool isMaxWeigth;
    GameObject newInternode;
    float annualLengthGrowth = Runge.runge(0.0f, 1.0f, 1.0f, .1f, new Runge.Function(Equation.equation)); //0.05f; // per ora lo fissiamo da codice, successivamente verra' calcolato
    float annualWidthGrowth = 0.01f;
    Vector3 annualGrowth;
    bool hasChild;


    // Start is called before the first frame update
    void Awake()
    {
        annualGrowth = new Vector3(annualWidthGrowth, annualLengthGrowth, annualWidthGrowth);
        setup = GameObject.Find("Setup").GetComponent<SetupScene>();
        timeManager = GameObject.Find("Time Manager").GetComponent<TimeManager>();
        currentLenght = 0f;
        currentWidth = 0f;
        geneticInfo = new Genetic(1);
    }

    // Update is called once per frame
    void Update()
    {
        Debug.Log("ciao, sono nato " + gameObject.name);
        Debug.Log(timeManager.IsYearElapsed + " " + gameObject.name);
        if (timeManager.IsYearElapsed)
        {

            //calcolo lunghezza e larghezza corrente
            float scaleFactorXZ = 1f;
            float scaleFactorY = 1f;
            if (gameObject.name != "Seed")
            {
                scaleFactorXZ = geneticInfo.InternodeMaxWidtht;
                scaleFactorY = geneticInfo.InternodeMaxLenght;
            }
            currentLenght = (gameObject.transform.lossyScale.y + annualGrowth.y);
            Debug.Log("Lunghezza corrente internodo = " + currentLenght);
            currentWidth = (gameObject.transform.lossyScale.x + annualGrowth.x);



            if(currentWidth > geneticInfo.InternodeMaxWidtht)
            {
                currentWidth = geneticInfo.InternodeMaxWidtht;
                isMaxWeigth = true;
            }

            if(currentLenght > geneticInfo.InternodeMaxLenght)
            {
                currentLenght = geneticInfo.InternodeMaxLenght;
                isMaxLenght = true;
            }



            gameObject.transform.localScale = new Vector3(currentWidth/scaleFactorXZ, currentLenght/scaleFactorY, currentWidth/scaleFactorXZ);



            if (isMaxLenght)
            {
                //crea nuovo internodo come figlio di quello corrente
                if (!hasChild && setup.Manager.InternodeNumber < geneticInfo.MaxInternodeNumber)
                {
                    setup.Manager.InternodeNumber++;
                    Debug.Log("nuovo internodo creato");
                    setup.LastInternodePosition = gameObject.transform.position + new Vector3(0f, geneticInfo.InternodeMaxLenght, 0f);
                    newInternode = Instantiate(setup.Prefab, new Vector3(setup.LastInternodePosition.x, setup.LastInternodePosition.y, setup.LastInternodePosition.z), Quaternion.identity);
                    newInternode.transform.localScale = Vector3.zero;
                    newInternode.transform.parent = setup.LastInternode.transform;
                    setup.LastInternode = newInternode;
                    //creare nuova classe gestore Gerarchia
                    newInternode.AddComponent<Growth>();
                    hasChild = true;
                }

                if (isMaxWeigth)
                {
                    gameObject.GetComponent<Growth>().enabled = false;
                }

            }
        }
    }

    public class Runge
    {
        //declare a delegate that takes a double and returns
        public delegate float Function(float t, float annualGrowth);

        public static float runge(float a, float b, float value, float step, Function f)
        {
            float t, w, k1, k2, k3, k4;
            t = a;
            w = value;
            for (int i = 0; i < (b - a) / step; i++)
            {
                k1 = step * f(t, w);
                k2 = step * f(t + step / 2f, w + k1 / 2f);
                k3 = step * f(t + step / 2f, w + k2 / 2f);
                k4 = step * f(t + step, w + k3);
                w = w + (k1 + 2f * k2 + 2f * k3 + k4) / 6f;
                t = a + i * step;
                //Console.WriteLine("{0} {1} ", Math.Round(t, 1), w);
            }

        }
    }

    class Equation
    {
        public static float equation(float t, float annualGrowth)
        {
            float y_prime;
            float k = 0.1f; //fattore di crescita, deve essere non hardcodato
            float maxHeight = 5.0f; //altezza max, deve essere non hardcodato
            return y_prime = k * annualGrowth * (1 - (annualGrowth / maxHeight));
        }

    }
}