Attributes[LSODA] = {HoldAllComplete, Protected}

LSODA[MethodData[v_, __]]["Version"] := v

LSODA[_]["StepMode"] := Automatic

LSODA[data_]["DenseOutputQ"] := DOQ[data]

LSODA /: InitializeMethod[LSODA, {stepmode_, denseoutput_}, nf_, ndstate_, rawmopts_] := Catch[Module[{denseQ, mdo, lsm, imps}, {lsm, mdo} = GetMethodOptionValues[LSODA, LSODAProperties, rawmopts]; CheckMaxDifferenceOrder[LSODA, mdo]; imps = TestImplicitSolver[LSODA, {"Newton", "LinearSolveMethod" -> lsm}]; denseQ = denseoutput === True; CreateMethodData[LSODA, {{}, False, mdo, imps, None, denseQ}]]]

Options[LSODA] = {"LinearSolveMethod" -> Automatic, "MaxDifferenceOrder" -> 12}