It sounds like wf1 and wf2 are supposed to be step inputs at t = 1 and t = 5. In Simulink these were implemented with Step blocks each followed by a Derivative block. Try replacing the Derivative blocks with Transfer Function blocks of the form:
s/(tau*s + 1)
where tau is suitably small. You might have to experiment with tau. As tau gets smaller, the solution will be more accurate but the simulation will run slower. Make sure the simulation is using a variable step solver and that the max step size is either auto or a small number. Might need to experiment with this parameter as well.
Another alternative would be to absorb the differentiation into downstream elements of the model. Wheter or not this is feasible depends on how the model is structured.
On the Matlab side, taking that numerical gradient might work. But the derivative of the step function is the Dirac delta function, so you can just as easily use the impulse function for each input, shift the outputs to the times the steps occur, and then sum the results.