function Propagator = createPropagator(x0,t0,propagator_data,tol,taylor_order)
    
    import java.io.*
    import java.lang.*
    import java.lang.Enum
    import java.text.*
    import java.util.*

    import org.orekit.attitudes.*
    import org.orekit.bodies.*
    import org.orekit.errors.OrekitException.*
    import org.orekit.forces.gravity.*
    import org.orekit.forces.gravity.potential.*
    import org.orekit.forces.gravity.potential.GravityFieldFactory.*
    import org.orekit.forces.maneuvers.*
    import org.orekit.frames.*
    import org.orekit.frames.LOFType.*
    import org.orekit.frames.LocalOrbitalFrame.*
    import org.orekit.frames.Frame.*
    import org.orekit.frames.FramesFactory
    import org.orekit.propagation.analytical.*
    import org.orekit.propagation.analytical.KeplerianPropagator
    import org.orekit.propagation.conversion.*
    import org.orekit.propagation.numerical.*
    import org.orekit.propagation.numerical.NumericalPropagator.*
    import org.orekit.propagation.semianalytical.dsst.*
    import org.orekit.propagation.semianalytical.dsst.forces.*
    import org.orekit.propagation.SpacecraftState
    import org.orekit.propagation.SpacecraftState.*
    import org.orekit.propagation.analytical.KeplerianPropagator.*
    import org.orekit.propagation.analytical.tle.*
    import org.orekit.orbits.*
    import org.orekit.orbits.KeplerianOrbit
    import org.orekit.orbits.KeplerianOrbit.*
    import org.orekit.orbits.Orbit.*
    import org.orekit.orbits.PositionAngle
    import org.orekit.orbits.PositionAngle.*
    import org.orekit.data.DataProvidersManager.*
    import org.orekit.data.ZipJarCrawler.*
    import org.orekit.time.*
    import org.orekit.time.AbsoluteDate.*
    import org.orekit.time.TimeScalesFactory.*
    import org.orekit.time.TimeScale.*
    import org.orekit.utils.*
    import org.hipparchus.*
    import org.hipparchus.geometry.euclidean.threed.*
    import org.orekit.models.earth.atmosphere.data.*
    import org.orekit.forces.gravity.potential.GravityFieldFactory.*
    
    crawler     = org.orekit.data.ZipJarCrawler('orekit-data.zip');
    % DM          = org.orekit.data.DataProvidersManager.getInstance(); % Not valid for OREKIT 10.3.x and after
    DM = org.orekit.data.DataContext.getDefault().getDataProvidersManager();
    DM.addProvider(crawler);
    
    % Initial data
    dx0     = x0;
    refDate = date_conversion(t0,'matlab2orekit');
    % endDate = date_conversion(end_time,'matlab2orekit');

    % Obtain constants
    mu  = Constants.WGS84_EARTH_MU;
    Re  = Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
    f   = Constants.WGS84_EARTH_FLATTENING;

    % Define constants for SRP
    dRef = 149597870000.0;
    pRef = 4.5605*1E-06;

    % Define frames and time
    ITRF    = FramesFactory.getITRF(ITRFVersion.ITRF_2000,IERSConventions.IERS_2010, 1); % ITRF_2000 as in the OEMs
    earth   = org.orekit.bodies.OneAxisEllipsoid(Re, f, ITRF);
    rotatingFrame = earth.getBodyFrame();
    inertialFrame = FramesFactory.getEME2000();
    UTC = TimeScalesFactory.getUTC();

    %%PERTURBATIONS FROM PROPAGATOR DATA
    % Prepare gravity model
    GravityFieldFactory.clearPotentialCoefficientsReaders();
    GravityFieldFactory.clearOceanTidesReaders();

    % GravityFieldFactory.addOceanTidesReader(FESCnmSnmReader('fes2004_Cnm-Snm.dat',1));
    GravityFieldFactory.addPotentialCoefficientsReader(ICGEMFormatReader('eigen-6s.gfc', true));
    gravity = GravityFieldFactory.getConstantNormalizedProvider(propagator_data.harmonics(1),propagator_data.harmonics(2));

    % isotropic model. This means nothing depends on attitude
    drag = org.orekit.forces.drag.IsotropicDrag(propagator_data.S,propagator_data.CD);

    if strcmp(propagator_data.drag,'Marshall')
        msafe = MarshallSolarActivityFutureEstimation(MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,javaMethod('valueOf','org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation$StrengthLevel', 'AVERAGE'));
        earthAtmosphere = org.orekit.models.earth.atmosphere.NRLMSISE00(msafe, org.orekit.bodies.CelestialBodyFactory.getSun(),earth);
        dragF = org.orekit.forces.drag.DragForce(earthAtmosphere,drag);

    elseif strcmp(propagator_data.drag,'HarrisPriester')
        earthAtmosphere = org.orekit.models.earth.atmosphere.HarrisPriester(org.orekit.bodies.CelestialBodyFactory.getSun(),earth);
        dragF = org.orekit.forces.drag.DragForce(earthAtmosphere,drag);

    else % strcmp(propagator_data.drag,'DTM2000')
        msafe = MarshallSolarActivityFutureEstimation(MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,javaMethod('valueOf','org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation$StrengthLevel', 'AVERAGE'));
        earthAtmosphere = org.orekit.models.earth.atmosphere.DTM2000(msafe, org.orekit.bodies.CelestialBodyFactory.getSun(),earth);
        dragF = org.orekit.forces.drag.DragForce(earthAtmosphere,drag);
    end

    % Complete perturbations with lunisolar and SRP if required
    sunF    = org.orekit.forces.gravity.ThirdBodyAttraction(org.orekit.bodies.CelestialBodyFactory.getSun());
    moonF   = org.orekit.forces.gravity.ThirdBodyAttraction(org.orekit.bodies.CelestialBodyFactory.getMoon());
    SRPCoef = 1.4;
    rad     = org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient(SRPCoef,propagator_data.SRPArea);
    SRPF    = org.orekit.forces.radiation.SolarRadiationPressure(dRef,pRef,org.orekit.bodies.CelestialBodyFactory.getSun(),Re,rad);

    % Integrator parameters
    minStep = tol.minStep;
    maxStep = tol.maxStep;
    initStep = tol.initStep;
    positionTolerance = tol.positionTolerance;

    % Set initial orbit
    if taylor_order >= 0
        params = 6;
        % Order of derivation of the DerivativeStructures
        order = taylor_order;
        factory = org.hipparchus.analysis.differentiation.DSFactory(params, order);
        x0 = factory.variable(0,dx0(1));
        y0 = factory.variable(1,dx0(2));
        z0 = factory.variable(2,dx0(3));
        vx0 = factory.variable(3,dx0(4));
        vy0 = factory.variable(4,dx0(5));
        vz0 = factory.variable(5,dx0(6));

        % Sometimes we will need the field of the DerivativeStructure to build new instances
        field = y0.getField();
        Fmu = factory.constant(mu);
        massF=factory.constant(propagator_data.mass);

        pos_0       = org.hipparchus.geometry.euclidean.threed.FieldVector3D(x0, y0, z0);
        vel_0       = org.hipparchus.geometry.euclidean.threed.FieldVector3D(vx0, vy0, vz0);
        pv_0        = org.orekit.utils.FieldPVCoordinates(pos_0, vel_0);
        FsimDate    = org.orekit.time.FieldAbsoluteDate(field, refDate);
        pvt_0       = org.orekit.utils.TimeStampedFieldPVCoordinates(FsimDate, pv_0);
        %Initialize the orbit
        orbit  = org.orekit.orbits.FieldCartesianOrbit(pvt_0, inertialFrame, Fmu);
    end
    pos_0NF       = org.hipparchus.geometry.euclidean.threed.Vector3D(dx0(1:3));
    vel_0NF       = org.hipparchus.geometry.euclidean.threed.Vector3D(dx0(4:6));
    pv_0NF        = org.orekit.utils.PVCoordinates(pos_0NF, vel_0NF);
    pvt_0NF       = org.orekit.utils.TimeStampedPVCoordinates(refDate, pv_0NF);
    %Initialize the orbit
    orbitNF       = org.orekit.orbits.CartesianOrbit(pvt_0NF, inertialFrame, mu);

    toleranc = org.orekit.propagation.numerical.NumericalPropagator.tolerances(positionTolerance, orbitNF, orbitNF.getType());
    absTolerance = toleranc(1,:);
    relTolerance = toleranc(2,:);

    if taylor_order >= 0
        Integrator = org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator(...
            field, minStep, maxStep, absTolerance, relTolerance);
    else
        Integrator = org.hipparchus.ode.nonstiff.DormandPrince853Integrator(...
            minStep, maxStep, absTolerance, relTolerance);
    end
    Integrator.setInitialStepSize(initStep);
%     Propagator
    if taylor_order >= 0
        Propagator = org.orekit.propagation.numerical.FieldNumericalPropagator(field, Integrator);
    else
        Propagator = org.orekit.propagation.numerical.NumericalPropagator(Integrator);
    end
    
    %   Initial State
    if taylor_order >= 0
        SS0 = org.orekit.propagation.FieldSpacecraftState(orbit,massF);
    else
        SS0 = org.orekit.propagation.SpacecraftState(orbitNF,propagator_data.mass);
    end
    Propagator.setInitialState(SS0);

    %   Atraction coefficient
    if taylor_order >= 0
        Propagator.setMu(Fmu);
    else
        Propagator.setMu(orbitNF.getMu());
    end
    %   Orbit Type
    Propagator.setOrbitType(orbitNF.getType());
    %   Add force models
    Propagator.addForceModel(org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel(rotatingFrame, gravity));
    if strcmp(propagator_data.drag,'none')
    else
        Propagator.addForceModel(dragF);
    end
    if propagator_data.sun==1
        Propagator.addForceModel(sunF);
    end
    if propagator_data.moon==1
        Propagator.addForceModel(moonF);
    end
    if propagator_data.SRP==1
        Propagator.addForceModel(SRPF);
    end

end