How do I implement the midpoint method when solving a pde?

I have one function that uses a for loop to iterate through every timestep to solve my pde and within that for loop I call three separate functions that solve the interdependent variables of my pde. I want to implement the midpoint method in one of these functions but I'm stuck because I'm not sure how to find the intermediate time step (t+h/2) - the whole function is being solved at one timestep.

I've tried the following method (simplified here)

itimestepvector = 1:dt:tf;
for itimestep = itimestepvector
    currentA = update_A_FE(state,...)   %update using forward euler
    lastB = currentB                    %I use the lastB in the midpoint method further down
    currentB = update_B_FE(state,...)

      %updating C using midpoint method
      k1 = dt*find_dcdt(state,currentc,currentB,...) %find_dcdt returns dC/dt using various constants and currentB
      k2 = dt*find_dcdt(state,currentc+0.5*dt*k1,0.5*(currentB+lastB))
      currentc = currentc + k2
  end

This method kind of works but is very unstable so presumably it's not actually the midpoint method. It also means C is one timestep behind B and A but I'm not too worried about that since my timestep is small and my tf very large.

Essentially, my problem is that find_dcdt only iterates through space not time (that's done in the overarching code above) and so I'm not sure how I'm supposed to use find_dcdt at t and t+dt/2 within only one timestep.