How to write the function of RK4 for 28*28 matrix size

I want to solve a first order ODE using RK4 method. My equation is dT/dt = C\(P(:,i)- G*(T-Tamb)). here C is diagonal capacitance matrix and size is [28*28] and G is a conductance matrix[28*28]. P is a power matrix [28*100] and I'm accessing each column in every time step. The part of my code is shown below. So I've three questions here.

1.Most of the solution I've found, first they solved the ode (homogeneous and particular) and they used these equations as function. In my case, do I need to solve my equation first or I can use it directly?

2. The value of K_1, K_2, K_3and K_4 are same in my case. But I'm expecting, they should be slightly different.

3. I tried the ODE45 solver. But the result is not as I expected. Any kind of suggestion is really helpful. Thank you in advance.

    Tamb = 75;                                              % Ambient temp
    h = 0.1;                                                % step size
    x = 1:h:10;
    y = zeros(28,length(x));                                % initialzation
    d0= [82.70,84,79.55,78.05,84.2,85.70,81.23,79.52,79.28,81.47,87.08,84.68,...
    77.73,79.47,84.24,82.79,82.79,82.79,82.79,82.79,82.79,82.79,82.79,82.79,82.79,82.79,82.79,82.79]; % initial value of y
    y(:,1) = d0;
    for i=1:(length(x)-1)                                   % calculation loop
        dydt = (C)\(P(:,i)-G*(y(:,i)-Tamb));                % Transient Equation
        F_xy = @(y)dydt ;
        k_1 = F_xy(y(:,i));                                 % calculating coefficient
        k_2 = F_xy(y(:,i)+0.5*h*k_1);
        k_3 = F_xy(y(:,i)+0.5*h*k_2);
        k_4 = F_xy(y(:,i)+k_3*h); 
        y(:,i+1) = y(:,i) + (h/6)*(k_1+2*k_2+2*k_3+k_4);  % main equation

    end

[JaeHoon Jang]

Hi, Md Shahidul Alam. I'm studying English. I hope your understanding.

I'm not found the C, P, and G matrix of your code.

So, Arbitrary matrices C, P, and G add your code.

    Tamb = 75;                                              % Ambient temp
    h = 0.1;                                                % step size
    x = 1:h:10;
    y = zeros(28,length(x));                                % initialzation
    d0= [82.70,84,79.55,78.05,84.2,85.70,81.23,79.52,79.28,81.47,87.08,84.68,...
    77.73,79.47,84.24,82.79,82.79,82.79,82.79,82.79,82.79,82.79,82.79,82.79,82.79,82.79,82.79,82.79]; % initial value of y
    y(:,1) = d0;

    V = ones(1,28);  %Arbitrary matrices
    C = diag(V);
    P = rand(28,100);
    G = C;*

    for i=1:(length(x)-1)                                   % calculation loop
        dydt = (C)\(P(:,i)-G*(y(:,i)-Tamb));                % Transient Equation
        F_xy = @(y)dydt ;
        k_1 = F_xy(y(:,i));                                 % calculating coefficient
        k_2 = F_xy(y(:,i)+0.5*h*k_1);
        k_3 = F_xy(y(:,i)+0.5*h*k_2);
        k_4 = F_xy(y(:,i)+k_3*h); 
        y(:,i+1) = y(:,i) + (h/6)*(k_1+2*k_2+2*k_3+k_4);  % main equation
    end

To find the values of K_1, K_2, K_3, and K_4.

[Md Shahidul Alam]

Hi JaeHoon Jang, Thank you for adding arbitrary C, P and G matrix. In my main matlab code I've these matrices. I've not added here because it'll be bulky. So, after adding these matrices now do you have the solution for my problems? Thank you.