Performance tips for differential equation RHS

It sucks to say, but when I tested it against the original code, I noticed the outputs were all wrong. Back to square one I think. How could I make the macro expansion work? I now tested with some code I actually understand but it’s just as slow as the code I began with (7s). Here is the new einsum function (which works)

function einsum_expansion!(coupling, coupling_matrix, coupling_terms, N)
    for j = 1:2
        for i = 1:N
            local s = zero(Float64)
            for r = 1:N
                @inbounds s += coupling_matrix[i, r] * coupling_terms[j, r]
            end
            @inbounds coupling[i, j] = s
        end
    end
    nothing
end